/[PAMELA software]/PamelaLevel2/src/PamLevel2.cpp
ViewVC logotype

Diff of /PamelaLevel2/src/PamLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Fri Jun 16 16:43:55 2006 UTC revision 1.37 by mocchiut, Wed Mar 28 15:29:14 2007 UTC
# Line 7  Line 7 
7   * Default constructor   * Default constructor
8   */   */
9  PamTrack::PamTrack(){  PamTrack::PamTrack(){
10      trk_track  = this->TrkTrack::GetTrkTrack();      trk_track  = 0;
11      calo_track = this->CaloTrkVar::GetCaloTrkVar();      calo_track = 0;
12      tof_track  = this->ToFTrkVar::GetToFTrkVar();      tof_track  = 0;
13        candeleteobj = 0;
14  };  };
15  //--------------------------------------  //--------------------------------------
16  //  //
# Line 19  PamTrack::PamTrack(){ Line 20  PamTrack::PamTrack(){
20   * Constructor   * Constructor
21   */   */
22  PamTrack::PamTrack(TrkTrack* t, CaloTrkVar* c, ToFTrkVar* o){  PamTrack::PamTrack(TrkTrack* t, CaloTrkVar* c, ToFTrkVar* o){
23      trk_track  = this->TrkTrack::GetTrkTrack();  
24      calo_track = this->CaloTrkVar::GetCaloTrkVar();      trk_track  = 0;
25      tof_track  = this->ToFTrkVar::GetToFTrkVar();      calo_track = 0;
26      if(t) *trk_track  = *t;      tof_track  = 0;
27      if(c) *calo_track = *c;  //     if(t)trk_track  = new TrkTrack(*t);
28      if(o) *tof_track  = *o;  //     if(c)calo_track = new CaloTrkVar(*c);
29    //     if(o)tof_track  = new ToFTrkVar(*o);
30        if(t)trk_track  = t;
31        if(c)calo_track = c;
32        if(o)tof_track  = o;
33        candeleteobj = 0;
34    
35  };  };
36    PamTrack::PamTrack(const PamTrack& track){
37    
38        TrkTrack   *t = track.trk_track;
39        CaloTrkVar *c = track.calo_track;
40        ToFTrkVar  *o = track.tof_track;
41    
42        trk_track  = 0;
43        calo_track = 0;
44        tof_track  = 0;
45        if(t)trk_track  = new TrkTrack(*t);
46        if(c)calo_track = new CaloTrkVar(*c);
47        if(o)tof_track  = new ToFTrkVar(*o);
48        candeleteobj = 1;
49        
50    }
51    void PamTrack::Clear(){
52    
53    //    cout << "PamTrack::Clear() "<<candeleteobj<<endl;
54        if(candeleteobj){
55            if(trk_track)  trk_track->TrkTrack::Clear();
56            if(calo_track) calo_track->CaloTrkVar::Clear();//???
57            if(tof_track)  tof_track->ToFTrkVar::Clear();//???
58        }else{
59            trk_track  = 0;
60            calo_track = 0;
61            tof_track  = 0;
62        }
63    }
64    void PamTrack::Delete(){
65    //    cout << "PamTrack::Delete() "<<candeleteobj<<endl;
66        if(candeleteobj){
67            if(trk_track)  {
68                trk_track->TrkTrack::Clear();
69                delete trk_track;
70            }
71            if(calo_track){
72                calo_track->CaloTrkVar::Clear();//???
73                delete calo_track;
74            }
75            if(tof_track){
76                tof_track->ToFTrkVar::Clear();//???
77                delete tof_track;
78            }
79        }else{
80            Clear();
81        }
82    }
83  //--------------------------------------  //--------------------------------------
84  //  //
85  //  //
86  //--------------------------------------  //--------------------------------------
87  /**  /**
88   * Constructor   * Default Constructor
89   */   */
90  PamLevel2::PamLevel2(){  PamLevel2::PamLevel2(){
91      trk_obj  = this->TrkLevel2::GetTrkLevel2();    Initialize();
92      calo_obj = this->CaloLevel2::GetCaloLevel2();  };
93      tof_obj  = this->ToFLevel2::GetToFLevel2();  
94      trig_obj = this->TrigLevel2::GetTrigLevel2();  /**
95      s4_obj   = this->S4Level2::GetS4Level2();   * Constructor with given dir, list and detectors
96      nd_obj   = this->NDLevel2::GetNDLevel2();   */
97      ac_obj   = this->AcLevel2::GetAcLevel2();  PamLevel2::PamLevel2(TString ddir,TString list,TString detlist){
98      Initialize();
99      GetPamTree(GetListOfLevel2Files(ddir,list),detlist);
100      GetRunTree(GetListOfLevel2Files(ddir,list));
101    };
102    
103    /**
104     * Constructor with given dir and list
105     */
106    PamLevel2::PamLevel2(TString ddir,TString list){
107      Initialize();
108      GetPamTree(GetListOfLevel2Files(ddir,list),"");
109      GetRunTree(GetListOfLevel2Files(ddir,list));
110    };
111    
112    
113    void PamLevel2::Initialize(){
114            
115    //     trk2_obj  = TrkLevel2::GetTrkLevel2();
116    //     trk1_obj  = TrkLevel1::GetTrkLevel1();
117    //     trkh_obj  = TrkHough::GetTrkHough();
118    //     calo1_obj = CaloLevel1::GetCaloLevel1();
119    //     calo2_obj = CaloLevel2::GetCaloLevel2();
120    //     tof_obj   = ToFLevel2::GetToFLevel2();
121    //     trig_obj  = TrigLevel2::GetTrigLevel2();
122    //     s4_obj    = S4Level2::GetS4Level2();
123    //     nd_obj    = NDLevel2::GetNDLevel2();
124    //     ac_obj    = AcLevel2::GetAcLevel2();
125    //     orb_obj   = OrbitalInfo::GetOrbitalInfo();
126        
127        h0_obj    = 0;
128        trk0_obj  = 0;
129        calo0_obj  = 0;
130    
131        trk2_obj  = 0;
132        trk1_obj  = 0;
133        trkh_obj  = 0;
134        calo1_obj = 0;
135        calo2_obj = 0;
136        tof_obj   = 0;
137        trig_obj  = 0;
138        s4_obj    = 0;
139        nd_obj    = 0;
140        ac_obj    = 0;
141        orb_obj   = 0;
142    
143        run_obj   = 0;//new GL_RUN();
144        soft_obj   = 0;// Emiliano
145        irun = -1LL;
146        runfirstentry = 0LL;
147        runlastentry = 0LL;
148    
149        l0_file = NULL;
150        l0_tree = NULL;
151        iroot   = -1;
152        dbc     = 0;
153    
154        prevshift = 0;
155        
156        run_tree = NULL;
157        run_tree_clone = NULL;
158        sel_tree = NULL;
159        sel_tree_clone = NULL;
160        
161        irunentry = -1LL;
162        pam_tree = NULL;
163        for(Int_t i=0; i<NCLONES; i++ )pam_tree_clone[i]=NULL;
164    
165        host = "mysql://localhost/pamelaprod";
166        user = "anonymous";
167        psw = "";
168        const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
169        const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
170        const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
171        if ( !pamdbhost ) pamdbhost = "";
172        if ( !pamdbuser ) pamdbuser = "";
173        if ( !pamdbpsw ) pamdbpsw = "";
174        if ( strcmp(pamdbhost,"") ) host = pamdbhost;
175        if ( strcmp(pamdbuser,"") ) user = pamdbuser;
176        if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
177    
178        
179    
180    //    sorted_tracks = 0;//new TRefArray();
181        
182        CAL0 = false;
183        CAL1 = true;
184        CAL2 = true;
185        TRK2 = true;
186        TRK1 = false;
187        TRKh = false;
188        TRKh = false;
189        TRG  = true;
190        TOF  = true;
191        TOF0 = false;
192        S4   = true;
193        ND   = true;
194        AC   = true;
195        ORB  = true;
196        
197        RUN  = true;
198    
199        SELLI = -1;
200    
201        tsorted=0;
202        timage=0;
203        
204        howtosort = "+CAL+TOF";
205        //howtosort = "+TOF";
206        sortthr = 100.;
207    
208    };
209    /**
210     * Delete the event (act as Dtor)
211     */
212    void PamLevel2::Delete(){
213            
214        if(run_obj)delete run_obj;
215        if(soft_obj)delete soft_obj; //Emiliano
216    
217    //    cout << "void PamLevel2::Clear()"<<endl;
218        if(h0_obj)   delete h0_obj;
219        if(trk0_obj) delete trk0_obj;      
220        if(calo0_obj) delete calo0_obj;    
221        if(trk1_obj) delete trk1_obj;      
222        if(trk2_obj) delete trk2_obj;
223        if(trkh_obj) delete trkh_obj;
224        if(calo1_obj)delete calo1_obj;
225        if(calo2_obj)delete calo2_obj;
226        if(tof_obj)  delete tof_obj;
227        if(trig_obj) delete trig_obj;
228        if(s4_obj)   delete s4_obj;
229        if(nd_obj)   delete nd_obj;
230        if(ac_obj)   delete ac_obj;
231        if(orb_obj)  delete orb_obj;
232        
233        if(tsorted){
234            tsorted->Delete();
235            delete tsorted;
236        }
237        if(timage){
238            timage->Delete();
239            delete timage;
240        }
241    
242        if(dbc){
243            dbc->Close();
244            delete dbc;
245        }
246        
247        if(l0_file)l0_file->Close();
248        //    if(pam_tree)pam_tree->Delete();;
249    
250      if ( pam_tree ){
251        //
252        // we have first to find which chains we have to delete, then delete the main chain and only after delete the friend chains. Any other order causes a segfault...
253        //
254        TList *temp = pam_tree->GetListOfFriends();
255        TList *contents  = new TList; // create chain friend list
256        contents->SetOwner();
257        TIter next(temp);
258        TChain *questo = 0;
259        while ( (questo = (TChain*) next()) ){
260          TString name =  questo->GetName();
261          contents->Add((TChain*)gROOT->FindObject(name.Data()));// add object to the list
262        };
263        //
264        // deleting the main chain
265        //
266        pam_tree->Delete();
267        //
268        // deleting the friends...
269        //
270        TIter next2(contents);
271        TChain *questa = 0;
272        while ( questa = (TChain*)next2() ){
273          TString name =  questa->GetName();
274          questa->Delete();      
275          questa=NULL;
276        };
277        //
278      };
279      pam_tree = NULL;
280    
281        if(run_tree)run_tree->Delete();;
282        if(sel_tree)sel_tree->Delete();;
283        for(Int_t i=0; i<NCLONES; i++ )if(pam_tree_clone[i])pam_tree_clone[i]->Delete();;
284        if(run_tree_clone)run_tree_clone->Delete();;
285        if(sel_tree_clone)sel_tree_clone->Delete();;
286        
287    
288    };
289    
290    /**
291     * Clear the event (NB! does not deallocate objects)
292     */
293    void PamLevel2::Clear(){
294            
295    //    cout << "void PamLevel2::Clear()"<<endl;
296    
297    //
298    // This method is called once for every entry but RunInfo and SoftInfo do not change until the next run so we cannot clear them here unless we don't
299    // want to load them for each event even if they are the same...
300    //
301    //    if(run_obj)delete run_obj;
302    //    if(run_obj) run_obj->Clear();  // Emiliano: Do not deallocate run_obj here, it will give segmentation fault! call clear instead
303    //    if(soft_obj) soft_obj->Clear();
304    
305        if(h0_obj)   h0_obj->Clear();      
306    //    if(trk0_obj) trk0_obj->Clear();  
307        if(calo0_obj) calo0_obj->Clear();  
308        if(trk1_obj) trk1_obj->Clear();    
309        if(trk2_obj) trk2_obj->Clear();
310        if(trkh_obj) trkh_obj->Clear();
311        if(calo1_obj)calo1_obj->Clear();
312        if(calo2_obj)calo2_obj->Clear();
313        if(tof_obj)  tof_obj->Clear();
314        if(trig_obj) trig_obj->Clear();
315        if(s4_obj)   s4_obj->Clear();
316        if(nd_obj)   nd_obj->Clear();
317        if(ac_obj)   ac_obj->Clear();
318        if(orb_obj)  orb_obj->Clear();
319        
320    //    if(sorted_tracks)sorted_tracks->Clear();
321    //    sorted_tracks.Clear();
322    
323        if(tsorted){
324            tsorted->Delete();
325        }
326        if(timage){
327            timage->Delete();
328        }
329    };
330    
331    void PamLevel2::Reset(){
332      //
333      // First of all clear everything
334      //
335      Clear();
336      //
337      // close and reset chains and pointers
338      //    
339      if ( pam_tree ){
340        //
341        // we have first to find which chains we have to delete, then delete the main chain and only after delete the friend chains. Any other order causes a segfault...
342        //
343        TList *temp = pam_tree->GetListOfFriends();
344        TList *contents  = new TList; // create chain friend list
345        contents->SetOwner();
346        TIter next(temp);
347        TChain *questo = 0;
348        while ( (questo = (TChain*) next()) ){
349          TString name =  questo->GetName();
350          contents->Add((TChain*)gROOT->FindObject(name.Data()));// add object to the list
351        };
352        //
353        // deleting the main chain
354        //
355        pam_tree->Delete();
356        //
357        // deleting the friends...
358        //
359        TIter next2(contents);
360        TChain *questa = 0;
361        while ( questa = (TChain*)next2() ){
362          TString name =  questa->GetName();
363          questa->Delete();      
364          questa=NULL;
365        };
366        //
367      };
368      pam_tree = NULL;
369      //
370      if(run_tree)run_tree->Delete();;
371      run_tree = NULL;
372      if(sel_tree)sel_tree->Delete();;
373      sel_tree = NULL;
374      //
375      // Close file
376      //
377      if(l0_file)l0_file->Close("R");
378      l0_file = NULL;
379      //
380      h0_obj    = 0;
381      trk0_obj  = 0;
382      calo0_obj  = 0;
383      //
384      trk2_obj  = 0;
385      trk1_obj  = 0;
386      trkh_obj  = 0;
387      calo1_obj = 0;
388      calo2_obj = 0;
389      tof_obj   = 0;
390      trig_obj  = 0;
391      s4_obj    = 0;
392      nd_obj    = 0;
393      ac_obj    = 0;
394      orb_obj   = 0;
395      //
396      // Reset run pointers
397      //
398      run_obj   = 0;//new GL_RUN();
399      soft_obj   = 0;// Emiliano
400      irun = -1;
401      runfirstentry = 0ULL;
402      runlastentry = 0ULL;    
403      //
404    };
405    
406    
407    //--------------------------------------
408    //
409    //
410    //--------------------------------------
411    void *PamLevel2::GetPointerTo(const char* c ){
412    
413        TString objname = c;
414    
415        if(!objname.CompareTo("TrkLevel1"))  {
416            if(!trk1_obj){
417                trk1_obj  = new TrkLevel1();
418                trk1_obj->Set();
419            }
420            return &trk1_obj;
421        };
422        if(!objname.CompareTo("TrkLevel2"))  {
423            if(!trk2_obj){  
424                trk2_obj  = new TrkLevel2();
425                trk2_obj->Set();
426            }
427            return &trk2_obj;
428        };
429        if(!objname.CompareTo("TrkHough"))   {
430            if(!trkh_obj)  trkh_obj  = new TrkHough();
431            return &trkh_obj;
432        };
433        if(!objname.CompareTo("CaloLevel1")) {
434            if(!calo1_obj) calo1_obj = new CaloLevel1();
435            return &calo1_obj;
436        };
437        if(!objname.CompareTo("CaloLevel2")) {
438            if(!calo2_obj){
439                calo2_obj = new CaloLevel2();
440                calo2_obj->Set();
441            };
442            return &calo2_obj;
443        };
444        if(!objname.CompareTo("ToFLevel2"))  {
445            if(!tof_obj){
446                tof_obj   = new ToFLevel2();
447                tof_obj->Set();
448            }
449            return &tof_obj;
450        };
451        if(!objname.CompareTo("TrigLevel2")) {
452            if(!trig_obj)  trig_obj  = new TrigLevel2();
453            return &trig_obj;
454        };
455        if(!objname.CompareTo("S4Level2"))   {
456            if(!s4_obj)    s4_obj    = new S4Level2();
457            return &s4_obj;
458        };
459        if(!objname.CompareTo("NDLevel2"))   {
460            if(!nd_obj)    nd_obj    = new NDLevel2();
461            return &nd_obj;
462        };
463        if(!objname.CompareTo("AcLevel2"))   {
464            if(!ac_obj)    ac_obj    = new AcLevel2();
465            return &ac_obj;
466        };
467        if(!objname.CompareTo("OrbitalInfo")){
468            if(!orb_obj)   orb_obj   = new OrbitalInfo();
469            return &orb_obj;
470        };
471        
472        if(!objname.CompareTo("RunInfo"))return &run_obj;
473    
474        if(!objname.CompareTo("SoftInfo"))return &soft_obj; // Emiliano
475    
476        return NULL;
477  };  };
478  //--------------------------------------  //--------------------------------------
479  //  //
480  //  //
481  //--------------------------------------  //--------------------------------------
482  /**  /**
483     * Retrieves the calorimeter track matching the seqno-th tracker stored track.
484     * (If seqno = -1 retrieves the self-trigger calorimeter track)
485     */
486     CaloTrkVar *PamLevel2::GetCaloStoredTrack(int seqno){
487        
488         if( !calo2_obj )return 0;
489    
490         if( calo2_obj->CaloLevel2::ntrk()==0 ){
491             cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no Calorimeter tracks are stored"<<endl;
492             return NULL;
493         };
494        
495         CaloTrkVar *c = 0;
496         Int_t it_calo=0;
497        
498         do{
499             c = calo2_obj->CaloLevel2::GetCaloTrkVar(it_calo);
500             it_calo++;
501         } while( c && seqno != c->trkseqno && it_calo < calo2_obj->CaloLevel2::ntrk());    
502        
503         if(!c || seqno != c->trkseqno){
504             c = 0;
505             if(seqno!=-1)cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match Calorimeter stored tracks"<<endl;
506         };
507         return c;
508        
509     };
510    //--------------------------------------
511     //
512     //
513    //--------------------------------------
514    /**
515      * Retrieves the ToF track matching the seqno-th tracker stored track.
516      * (If seqno = -1 retrieves the tracker-independent tof track)
517     */
518     ToFTrkVar *PamLevel2::GetToFStoredTrack(int seqno){
519            
520         if( !tof_obj )return 0;
521    
522         if( tof_obj->ToFLevel2::ntrk()==0 ){
523             cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no ToF tracks are stored"<<endl;
524             return NULL;
525         };
526        
527         ToFTrkVar *c = 0;
528         Int_t it_tof=0;
529        
530         do{
531             c = tof_obj->ToFLevel2::GetToFTrkVar(it_tof);
532             it_tof++;
533         } while( c && seqno != c->trkseqno && it_tof < tof_obj->ToFLevel2::ntrk());        
534        
535         if(!c || seqno != c->trkseqno){
536             c = 0;
537             if(seqno!=-1)cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match ToF stored tracks"<<endl;
538         };
539         return c;
540        
541     };
542    
543    //--------------------------------------
544     //
545     //
546    //--------------------------------------
547    /**
548     * Give the pamela track associated to a tracker track, retrieving related calorimeter and tof track information.
549     */
550     PamTrack* PamLevel2::GetPamTrackAlong(TrkTrack* t){
551        
552         cout <<"PamLevel2::GetPamTrackAlong(TrkTrack* t) **obsolete** "<<endl;
553         cout <<"(if you use it, remember to delete the PamTrack object)"<<endl;
554    
555         CaloTrkVar *c = 0;
556         ToFTrkVar  *o = 0;
557        
558         if(CAL2) c = GetCaloStoredTrack(t->GetSeqNo());
559         if(TOF) o = GetToFStoredTrack(t->GetSeqNo());
560        
561    //    if(t && c && o)track = new PamTrack(t,c,o);
562         PamTrack *track = new PamTrack(t,c,o);
563        
564         return track;
565    
566     };
567    //--------------------------------------
568    //
569    //
570    //--------------------------------------
571    /**
572   * Retrieves the it-th stored track.   * Retrieves the it-th stored track.
573   * It override TrkLevel2::GetTrack(int it).   * It override TrkLevel2::GetTrack(int it).
574   * @param itrk Track number, ranging from 0 to GetNTracks().   * @param itrk Track number, ranging from 0 to GetNTracks().
# Line 55  PamLevel2::PamLevel2(){ Line 576  PamLevel2::PamLevel2(){
576    
577  PamTrack* PamLevel2::GetStoredTrack(Int_t itrk){  PamTrack* PamLevel2::GetStoredTrack(Int_t itrk){
578            
579      // retrieve itrk-th tracker stored track      cout <<"PamLevel2::GetStoredTrack(Int_t itrk) **to-be-updated** "<<endl;
580      TrkTrack   *t = 0;      cout <<"for the moment, better use separately the methods: TrkLevel2::GetStoredTrack(seqno) CaloLevel2::GetCaloTrkVar(Int_t notrack) ToFLevel2::GetToFTrkVar(Int_t notrack)"<<endl;
581      CaloTrkVar *c = 0;      cout <<"(if you use it, remember to delete the PamTrack object)"<<endl;
582      ToFTrkVar  *o = 0;      PamTrack *track = 0;
583        
584      if( itrk >=0 && itrk < TrkLevel2::ntrk() ){      if( itrk >=0 && itrk < trk2_obj->TrkLevel2::ntrk() ){
585          t = TrkLevel2::GetStoredTrack(itrk);              
586          c = CaloLevel2::GetCaloTrkVar(t->GetSeqNo());          TrkTrack *t = trk2_obj->TrkLevel2::GetStoredTrack(itrk);
587          o = ToFLevel2::GetToFTrkVar(t->GetSeqNo());          track = GetPamTrackAlong(t);
588            
589            
590    
591        }else{
592            cout << "PamLevel2::GetStoredTrack(int) : tracker track SeqNo "<< itrk <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
593      };      };
     //  retrieve related ToF-track  
594            
     // hence create a "PamTrack" object  
     PamTrack *track = 0;  
     if(t && c && o)track = new PamTrack(t,c,o);  
595      return track;      return track;
596            
597  }  }
598  //--------------------------------------  //--------------------------------------
599  //  //
600    
601    /**
602     * Sort physical (tracker) tracks and stores them in the TRefArray (of TrkTrack objects) which pointer is  PamLevel2::sorted_tracks. Left here as backward compatibility method.
603     **/
604    void PamLevel2::SortTracks(TString how){
605      printf(" WARNING! obsolete, use SortTracks() and SetSortingMethod(TString) instead! \n Setting sorting method to %s \n",how.Data());
606      howtosort = how;  
607      SortTracks();
608    };
609    
610  //  //
611  //--------------------------------------  //--------------------------------------
612  /**  /**
613   * Sort physical (tracker) tracks and stores them in a TObjectArray (of TrkTrack objects).   * Sort physical (tracker) tracks and stores them in the TRefArray (of TrkTrack objects) which pointer is  PamLevel2::sorted_tracks.
614   * The total number of physical tracks is given by GetNTracks() and the it-th physical track can be retrieved by means of the method GetTrack(int it).   * @param how String to set the sorting cryterium (es: "CAL" or "TRK+CAL+TOF" ecc...).
615   * This method overrides TrkLevel2::GetTracks(), where sorting is done by decreasing number of fit points and increasing chi^2.   * Sorting cryteria:
616   * PamLevel2::GetTracks() keeps the same track order given by TrkLevel2::GetTracks(), but checks image selection by using calorimeter and ToF tracking information.   * TRK: lower chi**2
617     * CAL: lower Y spatial residual on the first calorimeter plane
618     * TOF: bigger numebr of hit PMTs along the track, on S12 S21 S32 (where paddles are along the Y axis).
619     * The default sorting cryterium is "TOF+CAL".
620     *
621     * The total number of physical tracks is always given by GetNTracks() and the it-th physical track can be retrieved by means of the methods GetTrack(int it) and GetTrack(int it, TString how).
622   */   */
623  TClonesArray *PamLevel2::GetTracks(){  void PamLevel2::SortTracks(){
624      TString how = howtosort;
625    
626      // define new array for sorted tracker tracks    //    cout <<" PamLevel2::SortTracks(TString how) "<<endl;
627      TClonesArray *aa = new TClonesArray("TrkTrack");    if( !trk2_obj ){
628      TClonesArray &sorted = *aa;      cout << "void PamLevel2::SortTracks():  TrkLevel2 not loaded !!!";
629        return;
630      // loop over the tracks sorted by the tracker    };
631      for(Int_t i=0; i < TrkLevel2::GetNTracks(); i++){    //    cout << "call SortTracs() "<<endl;
632      //Save current Object count
633          TrkTrack *ts=0;    Int_t ObjectNumber = TProcessID::GetObjectCount();
634    
635          // get tracker tracks    //    cout << "ObjectNumber  "<<ObjectNumber <<endl;
636          TrkTrack   *tp = TrkLevel2::GetTrack(i);                    //tracker  
637          CaloTrkVar *cp = CaloLevel2::GetCaloTrkVar(tp->GetSeqNo()); //calorimeter    //     if(!sorted_tracks)sorted_tracks = new TRefArray();
638          ToFTrkVar  *op = ToFLevel2::GetToFTrkVar(tp->GetSeqNo());   //tof    //     sorted_tracks->Clear();
639      //    sorted_tracks.Clear();
640          // if track has an image, check image selection  
641          if(tp->HasImage()){    if(!tsorted)tsorted = new TClonesArray("PamTrack",trk2_obj->GetNTracks());
642              TrkTrack   *ti = TrkLevel2::GetTrackImage(i);              //tracker (image)    tsorted->Delete();
643              CaloTrkVar *ci = CaloLevel2::GetCaloTrkVar(ti->GetSeqNo());//calorimeter (image)    TClonesArray &ttsorted = *tsorted;
644              ToFTrkVar  *oi = ToFLevel2::GetToFTrkVar(ti->GetSeqNo());  //tof (image)    if(!timage)timage = new TClonesArray("PamTrack",trk2_obj->GetNTracks());
645      timage->Delete();
646              //assign starting scores    TClonesArray &ttimage = *timage;
647              Int_t tp_score = 1;  
648              Int_t ti_score = 0;  
649              // ------------------------    // loop over the tracks sorted by the tracker
650              // calorimeter check    Bool_t use_TRK = how.Contains("TRK", TString::kIgnoreCase);
651              // ------------------------    Bool_t use_CAL = how.Contains("CAL", TString::kIgnoreCase);
652              if(    Bool_t use_TOF = how.Contains("TOF", TString::kIgnoreCase);
653                        npcfit[1] > 3   &&  //no. of fit planes on Y view      
654  //              varcfit[1] < 50.&&  //fit variance on Y view    if( !CAL2 &&  use_CAL) use_CAL = false;
655                        true){    if( !TOF &&  use_TOF) use_TOF = false;
656        
657                  Float_t resy_p = cp->tbar[0][1] - cbar[0][1];if(resy_p < 0)resy_p= - resy_p;    if( !TRK2 ){
658                  Float_t resy_i = ci->tbar[0][1] - cbar[0][1];if(resy_i < 0)resy_i= - resy_i;      //  cout << "SortTracks() : without tracker does not work!!! (not yet)" << endl;
659        return;
660                  if(resy_p <= resy_i) tp_score++;    };
                 else                 ti_score++;  
                       };  
             // ------------------------  
             // TOF check  
             // ------------------------      
   
             // ------------------------  
             // the winner is....  
             // ------------------------      
                       if(tp_score > ti_score) ts = tp;//the track sorted by the tracker!!  
                       else                    ts = ti;//its image!!  
661    
662          }else{    //   cout << "use_CAL "<<use_CAL<<" use_TOF "<<use_TOF<<" use_TRK "<<use_TRK <<endl;
663              ts = tp;      
664          };    for(Int_t i=0; i < trk2_obj->TrkLevel2::GetNTracks(); i++){
665                    
666          new(sorted[i]) TrkTrack(*ts); //save the track in the sorted array      TrkTrack   *ts = 0;
667      };      CaloTrkVar *cs = 0;
668        ToFTrkVar  *os = 0;
669            
670        // get tracker tracks
671        TrkTrack   *tp = trk2_obj->TrkLevel2::GetTrack(i);                    //tracker
672        CaloTrkVar *cp = GetCaloStoredTrack(tp->GetSeqNo());
673        ToFTrkVar  *op = GetToFStoredTrack(tp->GetSeqNo());
674    
675        TrkTrack   *ti = 0;              //tracker (image)
676        CaloTrkVar *ci = 0;
677        ToFTrkVar  *oi = 0;
678        //  cout << "trk track n. "<<i << " "<<hex<< tp <<dec<< endl;
679        // if track has an image, check image selection
680        if(tp->HasImage()){
681                
682          ti = trk2_obj->TrkLevel2::GetTrackImage(i);              //tracker (image)
683          ci = GetCaloStoredTrack(ti->GetSeqNo());
684          oi = GetToFStoredTrack(ti->GetSeqNo());
685                
686          //            cout << "its image "<<i << " "<<hex<< ti <<dec<< endl;
687    
688          //assign starting scores
689          Int_t tp_score = 0;  //main track sorted by the tracker
690          Int_t ti_score = 0;  //image track
691                
692          // ------------------------
693          // calorimeter check
694          // ------------------------
695          // check the Y spatial residual on the first calorimeter plane
696          // (cut on calorimeter variables from Emiliano)
697          if( use_CAL && !calo2_obj ){
698            cout << "void PamLevel2::SortTracks(): howtosort= "<<how<<" but CaloLevel2 not loaded !!!";
699            return;
700          };
701          if(
702             use_CAL            &&
703             calo2_obj->npcfit[1] > 5     &&   //no. of fit planes on Y view
704             calo2_obj->varcfit[1] < 1000. &&  //fit variance on Y view
705             cp && ci &&
706             true){
707    
708                    
709            Float_t resy_p = cp->tbar[0][1] - calo2_obj->cbar[0][1]; if(resy_p < 0)resy_p= - resy_p;
710            Float_t resy_i = ci->tbar[0][1] - calo2_obj->cbar[0][1]; if(resy_i < 0)resy_i= - resy_i;
711                    
712            if(resy_p <= resy_i) tp_score++;
713            else                 ti_score++;
714            
715            //      cout << "CALO "<<tp_score<<ti_score<<endl;
716    
717      return aa;        };
718          // ------------------------
719          // TOF check
720          // ------------------------          
721          // check the number of hit pmts along the track
722          // on S12 S21 and S32, where paddles are parallel to Y axis
723          if( use_TOF && !tof_obj ){
724            cout << "void PamLevel2::SortTracks(): howtosort= "<<how<<" but ToFLevel2 not loaded !!!";
725            return;
726          };
727          if( use_TOF && op && oi ){
728                    
729            //
730            Float_t sen = 0.;
731            for (Int_t ih=0; ih < op->npmtadc; ih++){
732              Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
733              if ( pl == 2 || pl == 3 || pl == 4 || pl == 5 ) sen += (op->dedx).At(ih);
734            };
735            for (Int_t ih=0; ih < oi->npmtadc; ih++){
736              Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
737              if ( pl == 2 || pl == 3 || pl == 4 || pl == 5 ) sen += (oi->dedx).At(ih);
738            };
739            //
740            if (  sen >= sortthr ){
741              //printf(" IS A NUCLEUS! en = %f \n",sen);
742              //
743              // is a nucleus use a different algorithm
744              //
745              Int_t nz = 6; Float_t zin[6];                                          // << define TOF z-coordinates
746              for(Int_t ip=0; ip<nz; ip++)
747                zin[ip] = tof_obj->ToFLevel2::GetZTOF(tof_obj->ToFLevel2::GetToFPlaneID(ip));     // << read ToF plane z-coordinates
748              Trajectory *tr = new Trajectory(nz,zin);
749              //
750              Int_t nphit_p =0;
751              Int_t nphit_i =0;
752              Float_t enhit_p = 0.;
753              Float_t enhit_i = 0.;
754              //
755              for (Int_t ih=0; ih < op->npmtadc; ih++){
756                Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
757                if(pl == 1 || pl == 2 || pl == 5){
758                  nphit_p++;
759                  enhit_p += (op->dedx).At(ih);
760                };
761              };
762              //
763              tp->DoTrack2(tr);
764              //
765              if ( fabs(tr->y[0]-oi->ytofpos[0]) < 2. ){
766                for (Int_t ih=0; ih < op->npmtadc; ih++){
767                  Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
768                  if(pl == 0){
769                    nphit_p++;
770                    enhit_p += (op->dedx).At(ih);
771                  };
772                };
773              };
774              if ( fabs(tr->y[3]-oi->ytofpos[1]) < 2. ){
775                for (Int_t ih=0; ih < op->npmtadc; ih++){
776                  Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
777                  if(pl == 3){
778                    nphit_p++;
779                    enhit_p += (op->dedx).At(ih);
780                  };
781                };
782              };
783              if ( fabs(tr->y[4]-oi->ytofpos[2]) < 2. ){
784                for (Int_t ih=0; ih < op->npmtadc; ih++){
785                  Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
786                  if(pl == 4){
787                    nphit_p++;
788                    enhit_p += (op->dedx).At(ih);
789                  };
790                };
791              };
792    
793              for (Int_t ih=0; ih < oi->npmtadc; ih++){
794                Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
795                if(pl == 1 || pl == 2 || pl == 5){
796                  nphit_i++;        
797                  enhit_i += (op->dedx).At(ih);
798                };
799              };
800              //
801              ti->DoTrack2(tr);
802              //
803              if ( fabs(tr->y[0]-oi->ytofpos[0]) < 2. ){
804                for (Int_t ih=0; ih < oi->npmtadc; ih++){
805                  Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
806                  if(pl == 0){
807                    nphit_i++;
808                    enhit_i += (op->dedx).At(ih);
809                  };
810                };
811              };
812              if ( fabs(tr->y[3]-oi->ytofpos[1]) < 2. ){
813                for (Int_t ih=0; ih < oi->npmtadc; ih++){
814                  Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
815                  if(pl == 3){
816                    nphit_i++;
817                    enhit_i += (op->dedx).At(ih);
818                  };
819                };
820              };
821              if ( fabs(tr->y[4]-oi->ytofpos[2]) < 2. ){
822                for (Int_t ih=0; ih < oi->npmtadc; ih++){
823                  Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
824                  if(pl == 4){
825                    nphit_i++;
826                    enhit_i += (op->dedx).At(ih);
827                  };
828                };
829              };            
830    
831                    
832              if(
833                 use_TOF            &&
834                 (nphit_p+nphit_i) !=0 &&  
835                 true){
836                        
837                //      printf(" seqno %i nphit_p %i nphit_i %i enhit_p %f enhit_i %f \n",trk2_obj->TrkLevel2::GetSeqNo(i),nphit_p,nphit_i,enhit_p,enhit_i);
838                //      printf(" score p %i score i %i \n",tp_score,ti_score);
839                //            if( enhit_p > enhit_i ) tp_score++;
840                //            if( nphit_p >= nphit_i && enhit_p > enhit_i ) tp_score++;
841                if ( nphit_p > nphit_i ) tp_score++;
842                if ( nphit_p < nphit_i ) ti_score++;
843                if ( nphit_p == nphit_i ){
844                  if ( enhit_p > enhit_i ) tp_score++;
845                  else ti_score++;
846                };
847                //      printf(" dopo score p %i score i %i \n",tp_score,ti_score);
848              };
849              delete tr;
850              //
851            } else {
852              //
853              // NOT a NUCLEUS
854              //
855              //printf(" NOT a NUCLEUS! en = %f \n",sen);
856    
857              Int_t nphit_p =0;
858              Int_t nphit_i =0;
859                    
860                    
861              /*                            cout << "track: npmtadc "<< op->npmtadc << endl;
862                cout << "track: npmttdc "<< op->npmttdc << endl;
863                cout << "image: npmtadc "<< oi->npmtadc << endl;
864                cout << "image: npmttdc "<< oi->npmttdc << endl;*/
865                    
866              for (Int_t ih=0; ih < op->npmtadc; ih++){
867                Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
868                if(pl == 1 || pl == 2 || pl == 5)nphit_p++;
869              };
870                    
871              for (Int_t ih=0; ih < oi->npmtadc; ih++){
872                Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
873                if(pl == 1 || pl == 2 || pl == 5)nphit_i++;
874              };
875                    
876              if(
877                 use_TOF            &&
878                 (nphit_p+nphit_i) !=0 &&  
879                 true){
880                        
881                if( nphit_p >= nphit_i) tp_score++;
882                else ti_score++;
883              };
884            };
885            //      cout << "TOF "<<tp_score<<ti_score<<endl;
886          };
887          if(tp_score == ti_score) use_TRK = true;
888          // ------------------------
889          // tracker check
890          // ------------------------
891          // chi**2 difference is not always large enough to distinguish among
892          // the real track and its image.
893          // Tracker check will be applied always when calorimeter and tof information is ambiguous.
894          if(use_TRK){
895            if(      tp->chi2 > 0 && tp->chi2 < ti->chi2 ) tp_score++ ;
896            else if( ti->chi2 > 0 && ti->chi2 < tp->chi2 ) ti_score++ ;
897            //      cout << "TRK "<<tp_score<<ti_score<<endl;
898          };
899                
900          // ------------------------
901          // the winner is....
902          // ------------------------          
903          if      (tp_score > ti_score) {
904            //              ts = tp;//the track sorted by the tracker!!
905            //              cs = cp;
906            //              os = op;
907            //      cout << "It's primary!" << endl;
908    
909          }else if (tp_score < ti_score) {
910    
911            //      cout << "It's image!" << endl;
912    
913            ts = ti;//its image!!
914            cs = ci;
915            os = oi;
916    
917            ti = tp;//its image!!
918            ci = cp;
919            oi = op;
920    
921            tp = ts;//its image!!
922            cp = cs;
923            op = os;
924    
925                    
926          }else {
927            //              ts = tp;
928            //              cs = cp;
929            //              os = op;
930            cout << "Warning - track image ambiguity not solved" << endl;
931            //                              cout << ts->GetNtot() << " "<< ts->chi2 << " " << npcfit[1] << " "<< nphit_p << endl;
932          };
933                
934        }else{
935          //            ts = tp;
936          //            cs = cp;
937          //            os = op;
938        };
939            
940        //  cout <<" SortTracks() "<<i<<" -- "<<ts<<endl;
941        //  sorted_tracks->Add(ts);//save the track in the sorted array
942        //  sorted_tracks.Add(ts);//save the track in the sorted array
943        //  sorted_tracks.Add(tp);//save the track in the sorted array
944        //  cout << "SortTracks:: sorted_tracks->Add(it) "<<i<<" "<<ts<<endl;
945        //  cout<<"o "<<tp<<endl;
946        //  cout<<"o "<<cp<<endl;
947        //  cout<<"o "<<op<<endl;
948        new(ttsorted[i]) PamTrack(tp,cp,op);
949        new(ttimage[i])  PamTrack(ti,ci,oi);
950      };
951    
952      if( tsorted->GetEntries() != trk2_obj->GetNTracks() ){
953        cout << "void PamLevel2::SortTracks(): tsorted->GetEntries() "<<tsorted->GetEntries()<<" != trk2_obj->GetNTracks() = "<<trk2_obj->GetNTracks() <<endl;
954        tsorted->Delete(); tsorted=0;
955        timage->Delete(); timage=0;
956      }
957    
958      //Restore Object count
959      //To save space in the table keeping track of all referenced objects
960      //We reset the object count to what it was at the beginning of the event.
961      TProcessID::SetObjectCount(ObjectNumber);
962        
963  };  };
964  //--------------------------------------  //--------------------------------------
965  //  //
966  //  //
967  //--------------------------------------  //--------------------------------------
968  /**  /**
969     * This method overrides TrkLevel2::GetTracks(), where sorting is done by decreasing number of fit points and increasing chi^2.
970     * PamLevel2::GetTracks() keeps the same track order given by TrkLevel2::GetTracks(), but checks image selection by using calorimeter and ToF tracking information.
971     */
972    // TRefArray *PamLevel2::GetTracks(){
973    
974    // //  *-*-*-*-*-*-*-*-*-*-*-*-*
975    //     SortTracks("+CAL+TOF");
976    // //  *-*-*-*-*-*-*-*-*-*-*-*-*
977    
978    // //   return  sorted_tracks;
979    //     return &sorted_tracks;
980        
981    //  };
982    TClonesArray *PamLevel2::GetTracks(){
983    
984    //  *-*-*-*-*-*-*-*-*-*-*-*-*
985        SortTracks();
986    //  *-*-*-*-*-*-*-*-*-*-*-*-*
987    
988        return tsorted;
989        
990     };
991    //--------------------------------------
992     //
993     //
994    //--------------------------------------
995    /**
996   * Retrieves the it-th Pamela "physical" track.   * Retrieves the it-th Pamela "physical" track.
997   * It override TrkLevel2::GetTrack(int it).   * It override TrkLevel2::GetTrack(int it).
998   * @param it Track number, ranging from 0 to GetNTracks().   * @param it Track number, ranging from 0 to GetNTracks().
999   */   */
1000  PamTrack *PamLevel2::GetTrack(int it){  PamTrack *PamLevel2::GetTrack(int it){
1001    
1002      TClonesArray *aa = this->GetTracks();  //    if(!trk2_obj) return 0;
1003      TClonesArray &sorted = *aa;          
1004    // //  *-*-*-*+-*-*-*-*-*-*-*-*-*
1005    //     SortTracks("+CAL+TOF");
1006    // //  *-*-*-*-*-*-*-*-*-*-*-*-*
1007    // //    if(!sorted_tracks)return 0;
1008    //     if(sorted_tracks.IsEmpty())return 0;
1009    
1010    //     PamTrack *track = 0;
1011        
1012    //     if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() && it<sorted_tracks.GetEntriesFast() ){    
1013    //      TrkTrack   *t = (TrkTrack*)sorted_tracks.At(it);
1014    //      track = GetPamTrackAlong(t);
1015    //     }else{
1016    //      cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1017    //     };
1018        
1019    //     return track;
1020    
1021      TrkTrack   *t = 0;  //    cout << "PamLevel2::GetTrack(int it) "<<endl;
1022      CaloTrkVar *c = 0;  //  *-*-*-*-*-*-*-*-*-*-*-*-*
1023      ToFTrkVar  *o = 0;      SortTracks();
1024    //  *-*-*-*-*-*-*-*-*-*-*-*-*
1025        if(!tsorted)return 0;
1026        if(!tsorted->GetEntries())return 0;
1027    
1028      if( it >=0 && it < TrkLevel2::GetNTracks() ){      PamTrack *track = 0;
1029          t = (TrkTrack*)sorted[it];      
1030          c = CaloLevel2::GetCaloTrkVar(t->GetSeqNo());      if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks()){
1031          o = ToFLevel2::GetToFTrkVar(t->GetSeqNo());  //      TrkTrack   *t = (TrkTrack*)sorted_tracks.At(it);
1032    //      track = GetPamTrackAlong(t);
1033            TClonesArray &t = *(tsorted);
1034            track = (PamTrack*)t[it];
1035            
1036        }else{
1037            cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1038      };      };
1039            
     // hence create a "PamTrack" object  
     PamTrack *track = 0;  
     if(t && c && o)track = new PamTrack(t,c,o);  
1040      return track;      return track;
1041        
1042  };  };
1043    
1044  //--------------------------------------  //--------------------------------------
1045  //  //
1046  //  //
1047  //--------------------------------------  //--------------------------------------
1048  /**  /**
1049   * Retrieves (if present) the image of the it-th Pamela "physical" track, sorted by the method PamLevel2::GetTracks().   * Retrieves (if present) the image of the it-th Pamela "physical" track, sorted by the method PamLevel2::SortTracks().
1050   * @param it Track number, ranging from 0 to GetNTracks().   * @param it Track number, ranging from 0 to GetNTracks().
1051   */   */
1052  PamTrack *PamLevel2::GetTrackImage(int it){  PamTrack *PamLevel2::GetTrackImage(int it){
1053    
1054      TClonesArray *aa = this->GetTracks();  //     if(!trk2_obj) return 0;
1055      TClonesArray &sorted = *aa;  
1056    // //  *-*-*-*-*-*-*-*-*-*-*-*-*
1057    //     SortTracks("+CAL+TOF");
1058    // //  *-*-*-*-*-*-*-*-*-*-*-*-*
1059    // //    if(!sorted_tracks)return 0;
1060    //     if(sorted_tracks.IsEmpty())return 0;
1061        
1062    //     PamTrack *image = 0;
1063        
1064    //     if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() ){
1065    //      TrkTrack *temp = (TrkTrack*)sorted_tracks.At(it);
1066    //      if( temp->HasImage() ){
1067    //          TrkTrack   *t = trk2_obj->TrkLevel2::GetStoredTrack(temp->GetImageSeqNo());
1068    //          image = GetPamTrackAlong(t);
1069    //      }else{
1070    //          cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
1071    //      };
1072    //     }else{
1073    //      cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1074    //     };
1075        
1076    //     return image;
1077    
     TrkTrack   *t = 0;  
     CaloTrkVar *c = 0;  
     ToFTrkVar  *o = 0;  
1078    
1079      if( it >=0 && it < TrkLevel2::GetNTracks() ){  //  *-*-*-*-*-*-*-*-*-*-*-*-*
1080          TrkTrack *temp = (TrkTrack*)sorted[it];      SortTracks();
1081          if( temp->HasImage() ){  //  *-*-*-*-*-*-*-*-*-*-*-*-*
1082              t = TrkLevel2::GetStoredTrack(temp->GetImageSeqNo());      if(!timage)return 0;
1083              c = CaloLevel2::GetCaloTrkVar(temp->GetImageSeqNo());      if(!timage->GetEntries())return 0;
1084              o = ToFLevel2::GetToFTrkVar(temp->GetImageSeqNo());      
1085        PamTrack *image = 0;
1086        
1087        if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() ){
1088            TClonesArray &t = *(tsorted);
1089            PamTrack *temp = (PamTrack*)t[it];
1090            if( temp->GetTrkTrack()->HasImage() ){
1091                TClonesArray &t = *(timage);
1092                image = (PamTrack*)t[it];
1093            }else{
1094    //          cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
1095          };          };
1096        }else{
1097            cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1098      };      };
1099            
     // hence create a "PamTrack" object  
     PamTrack *image = 0;  
     if(t && c && o)image = new PamTrack(t,c,o);  
1100      return image;      return image;
       
1101  }  }
1102    
1103  //--------------------------------------  //--------------------------------------
# Line 209  PamTrack *PamLevel2::GetTrackImage(int i Line 1105  PamTrack *PamLevel2::GetTrackImage(int i
1105  //  //
1106  //--------------------------------------  //--------------------------------------
1107  /**  /**
1108   * Get the Pamela detector trees and make them friends.   * Get the Pamela detector trees in a single file and make them friends.
1109     * @param f TFile pointer
1110     * @param detlist String to select trees to be included
1111     * @return Pointer to a TTree
1112   */   */
1113  TTree *PamLevel2::LoadPamTrees(TFile *f){  TTree *PamLevel2::GetPamTree(TFile *f, TString detlist ){
1114        
1115      TTree *Tout =0;  //     if( !detlist.IsNull() ) SetWhichTrees(detlist);
1116        //     else                    GetWhichTrees(f);
1117      // Tracker    if ( pam_tree ){
1118        printf("WARNING: TTree *PamLevel2::GetPamTree(TFile *fl, TString detlist) -- pam_tree already exists!\n ");
1119        return pam_tree;
1120      };
1121      //
1122    
1123        cout << "TTree *PamLevel2::GetPamTree(TFile *f, TString detlist ) -- obsolte "<<endl;
1124    
1125        SetWhichTrees(detlist);
1126    
1127        TTree *Trout =0;
1128    
1129        TString fname = f->GetName();
1130        if(!CheckLevel2File(fname))return NULL;
1131    
1132    //    UInt_t *found=0;    
1133    
1134        cout<< "GetPamTree(TFile*,TString): detector list --> ";
1135        if(TRK1)cout<<"TRK1 ";
1136        if(TRK2)cout<<"TRK2 ";
1137        if(TRKh)cout<<"TRKH ";
1138        if(CAL1)cout<<"CAL1 ";
1139        if(CAL2)cout<<"CAL2 ";
1140        if(TOF)cout<<"TOF ";
1141        if(TRG)cout<<"TRG ";
1142        if(AC)cout<<"AC ";
1143        if(ND)cout<<"ND ";
1144        if(S4)cout<<"S4 ";
1145        if(ORB)cout<<"ORB ";
1146        cout << endl;
1147        if(SELLI)cout<<">>> Found selection-list <<<"<<endl;
1148    
1149        f->cd();
1150    
1151    // Tracker
1152      TTree *T = (TTree*)f->Get("Tracker");      TTree *T = (TTree*)f->Get("Tracker");
1153      if(T) {      if(T && (TRK2||TRK1||TRKh)) {
1154          T->SetBranchAddress("TrkLevel2", GetPointerToTrk());          if(TRK2)T->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1155          cout << "Tracker      : set branch address TrkLevel2"<<endl;  //      else    T->SetBranchStatus("TrkLevel2",0,found);
1156          if(!Tout)Tout=T;          if(TRK2)cout << "Tracker      : set branch address TrkLevel2"<<endl;
1157            if(TRK1)T->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1158    //      else    T->SetBranchStatus("TrkLevel1",0,found);
1159            if(TRK1)cout << "Tracker      : set branch address TrkLevel1"<<endl;
1160            if(TRKh)T->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1161    //      else    T->SetBranchStatus("TrkHough",0,found);
1162            if(TRKh)cout << "Tracker      : set branch address TrkHough"<<endl;
1163            if(!Trout)Trout=T;
1164            else Trout->AddFriend(T);
1165      }else{      }else{
1166          cout << "Tracker      : missing tree"<<endl;          cout << "Tracker      : missing tree"<<endl;
1167      };      };
1168      // Calorimeter      // Calorimeter
1169      TTree *C = (TTree*)f->Get("Calorimeter");      TTree *C = (TTree*)f->Get("Calorimeter");
1170      if(C) {      if(C && (CAL2||CAL1)) {
1171          C->SetBranchAddress("CaloLevel2", GetPointerToCalo());          if(CAL2)C->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1172          cout << "Calorimeter  : set branch address CaloLevel2"<<endl;  //      else    C->SetBranchStatus("CaloLevel2",0,found);
1173          if(!Tout)Tout=C;          if(CAL2)cout << "Calorimeter  : set branch address CaloLevel2"<<endl;
1174          else Tout->AddFriend(C);          if(CAL1)C->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1175    //      else    C->SetBranchStatus("CaloLevel1",0,found);
1176            if(CAL1)cout << "Calorimeter  : set branch address CaloLevel1"<<endl;
1177            if(!Trout)Trout=C;
1178            else Trout->AddFriend(C);
1179      }else{      }else{
1180          cout << "Calorimeter  : missing tree"<<endl;          cout << "Calorimeter  : missing tree"<<endl;
1181      };      };
1182    
1183      // ToF          // ToF    
1184      TTree *O = (TTree*)f->Get("ToF");      TTree *O = (TTree*)f->Get("ToF");
1185      if(O) {      if(O && TOF) {
1186          O->SetBranchAddress("ToFLevel2", GetPointerToToF());          O->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1187          cout << "ToF          : set branch address ToFLevel2"<<endl;          cout << "ToF          : set branch address ToFLevel2"<<endl;
1188          if(!Tout)Tout=O;          if(!Trout)Trout=O;
1189          else Tout->AddFriend(O);          else Trout->AddFriend(O);
1190      }else{      }else{
1191          cout << "ToF         : missing tree"<<endl;          cout << "ToF         : missing tree"<<endl;
1192      };      };
1193      // Trigger      // Trigger
1194      TTree *R = (TTree*)f->Get("Trigger");      TTree *R = (TTree*)f->Get("Trigger");
1195      if(R) {      if(R && TRG) {
1196          R->SetBranchAddress("TrigLevel2", GetPointerToTrig());          R->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1197          cout << "Trigger      : set branch address TrigLevel2"<<endl;          cout << "Trigger      : set branch address TrigLevel2"<<endl;
1198          if(!Tout)Tout=O;          if(!Trout)Trout=R;
1199          else Tout->AddFriend(R);          else Trout->AddFriend(R);
1200      }else{      }else{
1201          cout << "Trigger      : missing tree"<<endl;          cout << "Trigger      : missing tree"<<endl;
1202      };      };
1203      // S4      // S4
1204      TTree *S = (TTree*)f->Get("S4");      TTree *S = (TTree*)f->Get("S4");
1205      if(S) {      if(S && S4) {
1206          S->SetBranchAddress("S4Level2", GetPointerToS4());          S->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1207          cout << "S4           : set branch address S4Level2"<<endl;          cout << "S4           : set branch address S4Level2"<<endl;
1208          if(!Tout)Tout=O;          if(!Trout)Trout=S;
1209          else Tout->AddFriend(S);          else Trout->AddFriend(S);
1210      }else{      }else{
1211          cout << "S4           : missing tree"<<endl;          cout << "S4           : missing tree"<<endl;
1212      };      };
1213      // Neutron Detector      // Neutron Detector
1214      TTree *N = (TTree*)f->Get("NeutronD");      TTree *N = (TTree*)f->Get("NeutronD");
1215      if(N) {      if(N && ND) {
1216          N->SetBranchAddress("NDLevel2", GetPointerToND());          N->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1217          cout << "NeutronD     : set branch address NDLevel2"<<endl;          cout << "NeutronD     : set branch address NDLevel2"<<endl;
1218          if(!Tout)Tout=O;          if(!Trout)Trout=N;
1219          else Tout->AddFriend(N);          else Trout->AddFriend(N);
1220      }else{      }else{
1221          cout << "NeutronD     : missing tree"<<endl;          cout << "NeutronD     : missing tree"<<endl;
1222      };      };
1223      // Anticounters      // Anticounters
1224      TTree *A = (TTree*)f->Get("Anticounter");      TTree *A = (TTree*)f->Get("Anticounter");
1225      if(A) {      if(A && AC) {
1226          A->SetBranchAddress("AcLevel2", GetPointerToAc());          A->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1227          cout << "Anticounter  : set branch address AcLevel2"<<endl;          cout << "Anticounter  : set branch address AcLevel2"<<endl;
1228          if(!Tout)Tout=O;          if(!Trout)Trout=A;
1229          else Tout->AddFriend(A);          else Trout->AddFriend(A);
1230      }else{      }else{
1231          cout << "Anticounter  : missing tree"<<endl;          cout << "Anticounter  : missing tree"<<endl;
1232      };      };
1233        // Orbital Info
1234        TTree *B = (TTree*)f->Get("OrbitalInfo");
1235        if(B && ORB) {
1236            B->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1237            cout << "OrbitalInfo  : set branch address OrbitalInfo"<<endl;
1238            if(!Trout)Trout=B;
1239            else Trout->AddFriend(B);
1240        }else{
1241            cout << "OrbitalInfo  : missing tree"<<endl;
1242        };
1243    
1244        TTree *L = (TTree*)f->Get("SelectionList");
1245        if(L && SELLI==1) {
1246    //      L->SetBranchAddress("RunEntry",&irun);
1247    //      cout << "SelectionList: set branch address RunEntry"<<endl;
1248    //      L->SetBranchAddress("EventEntry",&irunentry);
1249    //      cout << "SelectionList: set branch address EventEntry"<<endl;
1250    //      if(!Trout)Trout=O;
1251    //      else Trout->AddFriend("SelectionList");
1252            cout << " TTree *PamLevel2::GetPamTree(TFile, TString) >>> SelectionList not implemente!!!!"<<endl;
1253            sel_tree = 0;
1254        }else{
1255            cout << "SelectionList  : missing tree"<<endl;
1256        };
1257        
1258        cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
1259    
1260        pam_tree = (TChain*)Trout;
1261            
1262      return Tout;      return Trout;
1263            
1264  }  }
1265    //--------------------------------------
1266    //
1267    //
1268    //--------------------------------------
1269    /**
1270     * Get list of Level2 files.
1271     * @param ddir Level2 data directory.
1272     * @param flisttxt Name of txt file containing file list.
1273     * @return Pointer to a TList of TSystemFiles
1274     * If no input file list is given , all the Level2 files inside the directory are processed.
1275     */
1276    TList*  PamLevel2::GetListOfLevel2Files(TString ddir, TString flisttxt = ""){
1277            
1278        TString wdir = gSystem->WorkingDirectory();
1279        
1280        //    if(ddir=="")ddir = wdir;
1281    //      TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
1282        if ( ddir != ""){
1283          cout << "Level2 data directory : "<<  ddir << endl;
1284        } else {
1285          cout << "Level2 data directory not given as input: trying to evaluate from list or taking working directory " << endl;
1286        };
1287        TList *contents  = new TList; // create output list
1288        contents->SetOwner();
1289        
1290    //    char *fullpath;
1291    //    const char *fullpath;
1292        
1293        // if no input file list is given:  
1294        if ( flisttxt != "" ){
1295            
1296    //      if( !gSystem->IsFileInIncludePath(flisttxt,&fullpath) ){        
1297    //      if( !(fullpath = gSystem->FindFile("./",flisttxt)) ){
1298    //          cout <<"File "<<flisttxt<<" not found"<<endl;
1299    //          return 0;
1300    //      }      
1301    //      flisttxt = fullpath;
1302          if ( !flisttxt.EndsWith(".root") ){
1303    
1304            flisttxt = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
1305    
1306            if( !gSystem->ChangeDirectory(ddir) ){
1307                cout << "Cannot change directory : "<<ddir<<endl;
1308                return 0;
1309            }
1310            
1311            cout <<"Input file list : " << flisttxt <<endl;
1312            ifstream in;
1313            in.open(flisttxt, ios::in); //open input file list
1314            if(!in.good()){
1315                cout <<" ERROR opening the file "<<endl;
1316                gSystem->ChangeDirectory(wdir); // back to the working directory
1317                return 0;
1318            }      
1319            int line=0;
1320            while (1) {
1321                TString file;
1322                in >> file;
1323                if (!in.good()) break;
1324                line++;
1325    //          cout <<"(1) " << file << endl;
1326                if(file.IsNull()){
1327                    cout << "-- list interrupted at line "<<line <<endl;
1328                    break;
1329                }
1330                if(file.Contains("#"))file = file(0,file.First("#"));
1331    //          cout <<"(2) " << file << endl;
1332    //          if( gSystem->IsFileInIncludePath(file,&fullpath) ){
1333    //          if( (fullpath = gSystem->FindFile(ddir,file)) ){
1334                if( file.EndsWith(".root") ){
1335                  TString filedir;
1336                  if (ddir != ""){
1337                    filedir = ddir; // take the input dir
1338                  } else {
1339                    gSystem->ChangeDirectory(wdir); // back to the working directory
1340                    filedir = gSystem->DirName(file); // this will take the path if exist in the list otherwise it will return automatically the working dir
1341                  };
1342                    char *fullpath = gSystem->ConcatFileName(gSystem->DirName(filedir),gSystem->BaseName(file));
1343                    contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
1344                    delete fullpath;
1345                }
1346    //          }else{
1347    //              if(file.Data()!="")cout << "File: "<<file<<" ---> missing "<< endl;
1348    //          };
1349            };
1350            in.close();
1351          } else {
1352              if(flisttxt.Contains("#"))flisttxt = flisttxt(0,flisttxt.First("#"));
1353              char *fullpath = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
1354              contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
1355              delete fullpath;
1356          };        
1357        }else{
1358            
1359            cout << "No input file list given."<<endl;
1360            cout << "Check for existing root files."<<endl;
1361    //              cout << "Warking directory: "<< gSystem->WorkingDirectory()<< endl;
1362            if ( ddir == "" ){
1363              ddir = wdir;
1364              cout << "Level2 data directory : "<<  ddir << endl;
1365            };
1366    
1367            TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
1368            TList *temp = datadir->GetListOfFiles();
1369    //              temp->Print();
1370    //              cout << "*************************************" << endl;
1371            
1372            TIter next(temp);
1373            TSystemFile *questo = 0;
1374            
1375            
1376            while ( (questo = (TSystemFile*) next()) ) {
1377                TString name =  questo-> GetName();
1378                if( name.EndsWith(".root") ){
1379    //              const char *fullpath = gSystem->FindFile(ddir,name);
1380    //              char *fullpath;
1381    //              gSystem->IsFileInIncludePath(name,&fullpath);
1382                    char *fullpath = gSystem->ConcatFileName(gSystem->DirName(ddir),gSystem->BaseName(name));
1383                    contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));
1384                    delete fullpath;
1385                };
1386            }
1387            delete temp;
1388            delete datadir;
1389            
1390        };
1391        gSystem->ChangeDirectory(wdir); // back to the working directory
1392    //      cout << endl << "Selected files:" << endl;
1393    //      contents->Print();
1394        cout << contents->GetEntries()<<" files \n";
1395    //      cout << endl;
1396    //      cout << "Working directory: "<< gSystem->WorkingDirectory()<< endl;
1397        return contents;
1398    };
1399    //--------------------------------------
1400    //
1401    //
1402    //--------------------------------------
1403    /**
1404     * Get the Pamela detector chains from a list of files and make them friends.
1405     * @param fl Pointer to a TList of TSystemFiles
1406     * @param detlist String to select trees to be included
1407     * @return Pointer to a TChain
1408     */
1409    TChain *PamLevel2::GetPamTree(TList *fl, TString detlist ){
1410      //
1411      //
1412      //
1413      if ( pam_tree ){
1414        printf("WARNING: TChain *PamLevel2::GetPamTree(TList *fl, TString detlist) -- pam_tree already exists!\n ");
1415        return pam_tree;
1416      };
1417      //
1418        
1419        TChain *Trout =0;
1420    
1421    //    if( !detlist.IsNull() )SetWhichTrees(detlist);
1422        SetWhichTrees(detlist);
1423        
1424        cout<< "GetPamTree(TList*,TString): input detector list --> ";
1425        if(TRK1)cout<<"TRK1 ";
1426        if(TRK2)cout<<"TRK2 ";
1427        if(TRKh)cout<<"TRKH ";
1428        if(CAL1)cout<<"CAL1 ";
1429        if(CAL2)cout<<"CAL2 ";
1430        if(TOF)cout<<"TOF ";
1431        if(TRG)cout<<"TRG ";
1432        if(AC)cout<<"AC ";
1433        if(ND)cout<<"ND ";
1434        if(S4)cout<<"S4 ";
1435        if(ORB)cout<<"ORB ";
1436        cout << endl;
1437    
1438         TChain *T = 0;    
1439         TChain *C = 0;
1440         TChain *O = 0;
1441         TChain *R = 0;
1442         TChain *S = 0;
1443         TChain *N = 0;
1444         TChain *A = 0;
1445         TChain *B = 0;
1446    
1447         TChain *L = 0;
1448        
1449        if(TRK2||TRK1||TRKh) T = new TChain("Tracker");    
1450        if(CAL2||CAL1)       C = new TChain("Calorimeter");
1451        if(TOF)              O = new TChain("ToF");
1452        if(TRG)              R = new TChain("Trigger");
1453        if(S4)               S = new TChain("S4");
1454        if(ND)               N = new TChain("NeutronD");
1455        if(AC)               A = new TChain("Anticounter");
1456        if(ORB)              B = new TChain("OrbitalInfo");
1457         L = new TChain("SelectionList");
1458        
1459        // loop over files and create chains        
1460        TIter next(fl);
1461        TSystemFile *questo = 0;
1462        while ( (questo = (TSystemFile*) next()) ) {
1463            TString name =  questo->GetName();
1464            cout << "File: "<< name << endl;
1465            if( CheckLevel2File(name) ){
1466                if(TRK2||TRK1||TRKh) T->Add(name);
1467                if(CAL1||CAL2)       C->Add(name);
1468                if(TOF)              O->Add(name);
1469                if(TRG)              R->Add(name);
1470                if(S4)               S->Add(name);
1471                if(ND)               N->Add(name);
1472                if(AC)               A->Add(name);
1473                if(ORB)              B->Add(name);
1474                if(SELLI==1)         L->Add(name);
1475            };
1476        };
1477        
1478        cout << "done chain \n";
1479        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
1480    
1481    //    UInt_t *found=0;
1482    // Tracker
1483        if(T && (TRK2||TRK1||TRKh)) {
1484            if(TRK2)T->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1485    //      else    T->SetBranchStatus("TrkLevel2",0,found);
1486            if(TRK2)cout << "Tracker      : set branch address TrkLevel2"<<endl;
1487            if(TRK1)T->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1488    //      else    T->SetBranchStatus("TrkLevel1",0,found);
1489            if(TRK1)cout << "Tracker      : set branch address TrkLevel1"<<endl;
1490            if(TRKh)T->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1491    //      else    T->SetBranchStatus("TrkHough",0,found);
1492            if(TRKh)cout << "Tracker      : set branch address TrkHough"<<endl;
1493            if(!Trout)Trout=T;
1494            else Trout->AddFriend("Tracker");
1495        }else{
1496            cout << "Tracker      : missing tree"<<endl;
1497        };
1498        // Calorimeter
1499        if(C && (CAL2||CAL1)) {
1500            if(CAL2)C->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1501    //      else    C->SetBranchStatus("CaloLevel2",0,found);
1502            if(CAL2)cout << "Calorimeter  : set branch address CaloLevel2"<<endl;
1503            if(CAL1)C->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1504    //      else    C->SetBranchStatus("CaloLevel1",0,found);
1505            if(CAL1)cout << "Calorimeter  : set branch address CaloLevel1"<<endl;
1506            if(!Trout)Trout=C;
1507            else Trout->AddFriend("Calorimeter");
1508        }else{
1509            cout << "Calorimeter  : missing tree"<<endl;
1510        };
1511        // ToF    
1512        if(O && TOF) {
1513            O->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1514            cout << "ToF          : set branch address ToFLevel2"<<endl;
1515            if(!Trout)Trout=O;
1516            else Trout->AddFriend("ToF");
1517        }else{
1518            cout << "ToF         : missing tree"<<endl;
1519        };
1520        // Trigger
1521        if(R && TRG) {
1522            R->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1523            cout << "Trigger      : set branch address TrigLevel2"<<endl;
1524            if(!Trout)Trout=O;
1525            else Trout->AddFriend("Trigger");
1526        }else{
1527            cout << "Trigger      : missing tree"<<endl;
1528        };
1529        // S4
1530        if(S && S4) {
1531            S->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1532            cout << "S4           : set branch address S4Level2"<<endl;
1533            if(!Trout)Trout=O;
1534            else Trout->AddFriend("S4");
1535        }else{
1536            cout << "S4           : missing tree"<<endl;
1537        };
1538        // Neutron Detector
1539        if(N && ND) {
1540            N->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1541            cout << "NeutronD     : set branch address NDLevel2"<<endl;
1542            if(!Trout)Trout=O;
1543            else Trout->AddFriend("NeutronD");
1544        }else{
1545            cout << "NeutronD     : missing tree"<<endl;
1546        };
1547        // Anticounters
1548        if(A && AC) {
1549            A->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1550            cout << "Anticounter  : set branch address AcLevel2"<<endl;
1551            if(!Trout)Trout=O;
1552            else Trout->AddFriend("Anticounter");
1553        }else{
1554            cout << "Anticounter  : missing tree"<<endl;
1555        };
1556        // Orbital Info
1557        if(B && ORB) {
1558            B->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1559            cout << "OrbitalInfo  : set branch address OrbitalInfo"<<endl;
1560            if(!Trout)Trout=O;
1561            else Trout->AddFriend("OrbitalInfo");
1562        }else{
1563            cout << "OrbitalInfo  : missing tree"<<endl;
1564        };
1565        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
1566    
1567        // Selection List
1568        if(L && SELLI==1) {
1569            cout<<">>> Found selection-list <<<"<<endl;
1570            L->SetBranchAddress("RunEntry",&irun);
1571            cout << "SelectionList: set branch address RunEntry"<<endl;
1572            L->SetBranchAddress("EventEntry",&irunentry);
1573            cout << "SelectionList: set branch address EventEntry"<<endl;
1574            sel_tree = L;
1575    //      if(!Trout)Trout=O;
1576    //      else Trout->AddFriend("SelectionList");
1577            cout << "----------------------------------------------------" <<endl;
1578        }else{
1579    //      cout << "SelectionList  : missing tree"<<endl;
1580            if(L)L->Delete();
1581        };
1582        
1583    //    cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
1584    
1585        pam_tree = Trout;
1586    
1587        return Trout;
1588    }
1589    
1590    
1591    
1592    //--------------------------------------
1593    //
1594    //
1595    //--------------------------------------
1596    /**
1597     * Set branch addresses for Pamela friend trees
1598     */
1599    void PamLevel2::SetBranchAddress(TTree *t){
1600    
1601        TRK2    = TRK2 & t->GetBranchStatus("TrkLevel2");
1602        TRK1    = TRK1 & t->GetBranchStatus("TrkLevel1");
1603        TRKh    = TRKh & t->GetBranchStatus("TrkHough");
1604        CAL2    = CAL2 & t->GetBranchStatus("CaloLevel2");
1605        CAL1    = CAL1 & t->GetBranchStatus("CaloLevel1");
1606        TOF    = TOF & t->GetBranchStatus("ToFLevel2");
1607        TRG    = TRG & t->GetBranchStatus("TrigLevel2");
1608        S4     = S4  & t->GetBranchStatus("S4Level2");
1609        ND     = ND  & t->GetBranchStatus("NDLevel2");
1610        AC     = AC  & t->GetBranchStatus("AcLevel2");
1611        ORB    = ORB & t->GetBranchStatus("OrbitalInfo");
1612    
1613    
1614        // Tracker
1615        if(TRK1) {
1616            t->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel2"));
1617            cout << "Tracker      : set branch address TrkLevel1"<<endl;
1618        };
1619        if(TRK2) {
1620            t->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel1"));
1621            cout << "Tracker      : set branch address TrkLevel2"<<endl;
1622        };
1623        if(TRKh) {
1624            t->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1625            cout << "Tracker      : set branch address TrkHough"<<endl;
1626        };
1627        
1628        // Calorimeter
1629        if(CAL1) {
1630            t->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1631            cout << "Calorimeter  : set branch address CaloLevel1"<<endl;
1632        };
1633        if(CAL2) {
1634            t->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1635            cout << "Calorimeter  : set branch address CaloLevel2"<<endl;
1636        };
1637        
1638        // ToF    
1639        if(TOF) {
1640            t->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1641            cout << "ToF          : set branch address ToFLevel2"<<endl;
1642        };
1643        // Trigger
1644        if(TRG) {
1645            t->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1646            cout << "Trigger      : set branch address TrigLevel2"<<endl;
1647        };
1648        // S4
1649        if(S4) {
1650            t->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1651            cout << "S4           : set branch address S4Level2"<<endl;
1652        };
1653        // Neutron Detector
1654        if(ND) {
1655            t->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1656            cout << "NeutronD     : set branch address NDLevel2"<<endl;
1657        };
1658        // Anticounters
1659        if(AC) {
1660            t->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1661            cout << "Anticounter  : set branch address AcLevel2"<<endl;
1662        };
1663        // OrbitalInfo
1664        if(ORB) {
1665            t->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1666            cout << "OrbitalInfo  : set branch address OrbitalInfo"<<endl;
1667        };
1668        // SelectionList
1669        if(SELLI==1) {
1670            t->SetBranchAddress("RunEntry", &irun);
1671            cout << "SelectionList: set branch address RunEntry"<<endl;
1672            t->SetBranchAddress("EventEntry", &irunentry);
1673            cout << "SelectionList: set branch address EventEntry"<<endl;
1674        };
1675        
1676    }
1677    /**
1678     * Set branch addresses for Pamela friend trees
1679     */
1680    void PamLevel2::SetBranchAddress(TChain *t){
1681    
1682        TRK2    = TRK2 & t->GetBranchStatus("TrkLevel2");
1683        TRK1    = TRK1 & t->GetBranchStatus("TrkLevel1");
1684        TRKh    = TRKh & t->GetBranchStatus("TrkHough");
1685        CAL1    = CAL1 & t->GetBranchStatus("CaloLevel1");
1686        CAL2    = CAL2 & t->GetBranchStatus("CaloLevel2");
1687        TOF    = TOF & t->GetBranchStatus("ToFLevel2");
1688        TRG    = TRG & t->GetBranchStatus("TrigLevel2");
1689        S4     = S4  & t->GetBranchStatus("S4Level2");
1690        ND     = ND  & t->GetBranchStatus("NDLevel2");
1691        AC     = AC  & t->GetBranchStatus("AcLevel2");
1692        ORB    = ORB & t->GetBranchStatus("OrbitalInfo");
1693    
1694        // Tracker
1695         if(TRK2) {
1696            t->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1697            cout << "Tracker      : set branch address TrkLevel2"<<endl;
1698        };
1699        if(TRK1) {
1700            t->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1701            cout << "Tracker      : set branch address TrkLevel1"<<endl;
1702        };
1703         if(TRKh) {
1704            t->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1705            cout << "Tracker      : set branch address TrkHough"<<endl;
1706         };
1707        
1708        // Calorimeter
1709        if(CAL2) {
1710            t->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1711            cout << "Calorimeter  : set branch address CaloLevel2"<<endl;
1712        };  
1713        if(CAL1) {
1714            t->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1715            cout << "Calorimeter  : set branch address CaloLevel1"<<endl;
1716        };
1717        
1718        // ToF    
1719        if(TOF) {
1720            t->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1721            cout << "ToF          : set branch address ToFLevel2"<<endl;
1722        };
1723        // Trigger
1724        if(TRG) {
1725            t->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1726            cout << "Trigger      : set branch address TrigLevel2"<<endl;
1727        };
1728        // S4
1729        if(S4) {
1730            t->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1731            cout << "S4           : set branch address S4Level2"<<endl;
1732        };
1733        // Neutron Detector
1734        if(ND) {
1735            t->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1736            cout << "NeutronD     : set branch address NDLevel2"<<endl;
1737        };
1738        // Anticounters
1739        if(AC) {
1740            t->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1741            cout << "Anticounter  : set branch address AcLevel2"<<endl;
1742        };
1743        // OrbitalInfo
1744        if(ORB) {
1745            t->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1746            cout << "OrbitalInfo  : set branch address OrbitalInfo"<<endl;
1747        };
1748        // SelectionList
1749        if(SELLI==1) {
1750            t->SetBranchAddress("RunEntry", &irun);
1751            cout << "SelectionList: set branch address RunEntry"<<endl;
1752            t->SetBranchAddress("EventEntry", &irunentry);
1753            cout << "SelectionList: set branch address EventEntry"<<endl;
1754        };
1755        
1756    }
1757    
1758    
1759    //--------------------------------------
1760    //
1761    //
1762    //--------------------------------------
1763    /**
1764     * Get the Run tree chain from a list of files.
1765     * @param fl Pointer to a TList of TSystemFiles
1766     * @return Pointer to a TChain
1767     */
1768    TChain *PamLevel2::GetRunTree(TList *fl){
1769      //
1770      //
1771      //
1772      if ( run_tree ){
1773        printf("WARNING: TChain *PamLevel2::GetRunTree(TList *fl) -- run_tree already exists!\n ");
1774        return run_tree;
1775      };    
1776      //
1777      TChain *R = new TChain("Run");
1778        
1779        // loop over files and create chains        
1780        TIter next(fl);
1781        TSystemFile *questo = 0;
1782        while ( (questo = (TSystemFile*) next()) ) {
1783            TString name =  questo->GetName();
1784    //              cout << "File: "<< name << endl;
1785            if( CheckLevel2File(name) ){
1786                R->Add(name);
1787            };
1788        }
1789    
1790        if(R->GetNtrees()){
1791        
1792            R->SetBranchAddress("RunInfo", GetPointerTo("RunInfo"));
1793            cout << "Run         : set branch address RunInfo"<<endl;
1794            R->SetBranchAddress("SoftInfo", GetPointerTo("SoftInfo")); // Emiliano
1795            cout << "Software    : set branch address SoftInfo"<<endl; // Emiliano
1796        }else{
1797            delete R;
1798            R=0;
1799        }
1800    
1801        run_tree = R;
1802    
1803        return R;
1804        
1805    }
1806    //--------------------------------------
1807    //
1808    //
1809    //--------------------------------------
1810    /**
1811     * Get the Run tree  from a file.
1812     * @param f Pointer to a TFile
1813     * @return Pointer to a TTree
1814     */
1815    TTree *PamLevel2::GetRunTree(TFile *f){
1816      if ( run_tree ){
1817        printf("WARNING: TTree *PamLevel2::GetRunTree(TFile *f) -- run_tree already exists!\n ");
1818        return run_tree;
1819      };
1820            
1821    
1822        cout << "TTree *PamLevel2::GetRunTree(TFile *f) -- obsolte "<<endl;
1823    
1824        TTree *T = (TTree*)f->Get("Run");
1825        
1826        if(T){
1827            T->SetBranchAddress("RunInfo", GetPointerTo("RunInfo"));
1828            cout << "Run         : set branch address RunInfo"<<endl;
1829            T->SetBranchAddress("SoftInfo", GetPointerTo("SoftInfo")); // Emiliano
1830            cout << "Software    : set branch address SoftInfo"<<endl; // Emiliano
1831        }
1832    
1833        run_tree = (TChain*)T;
1834    
1835        return T;
1836        
1837    }
1838    /**
1839     * Update the runinfo informations (to be used to have Run infos event by event basis)
1840     * @param run Pointer to the chain/tree which contains run infos
1841     * @return true if a new run has been read, false if it is still the same run
1842     */
1843    Bool_t PamLevel2::UpdateRunInfo(TChain *run, Long64_t iev){
1844      //
1845      // check if we have already called once GetEntry, if not call it
1846      //
1847        cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, Long64_t iev) --- NON FUNZIONA BENISSIMO.... "<<endl;
1848        if(!run){
1849              cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing RunInfo tree "<<endl;
1850              return(false);        
1851        }
1852        if ( run->GetEntries() <= 0 ) return(false);
1853      //
1854      
1855      //  Int_t oldrun = irun;
1856      Long64_t oldrun = irun;
1857    
1858      // --------------------------------------
1859      // if it is a full file (not preselected)
1860      // --------------------------------------
1861      if(SELLI==0){
1862    
1863          //
1864          // the absolute time is necessary to relate the event with the run
1865          //
1866          if( !GetOrbitalInfo() ){
1867              cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing OrbitalInfo "<<endl;
1868              return(false);
1869          }
1870    
1871          //
1872          // the first time the routine is called, set run search from the beginning
1873          //
1874          if ( irun < 0LL ){
1875              irun = 0LL;
1876              run->GetEntry(irun);
1877              runfirstentry = 0LL;
1878              runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1879              if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runlastentry -= 1LL;
1880          };      
1881          //
1882           if ( irun == run->GetEntries()-1LL &&
1883               !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME &&
1884                 GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)
1885                ){
1886             irun = -1LL;
1887             runfirstentry = 0LL;
1888             runlastentry = -1LL;
1889           };
1890          // modificato il controllo sull'aggiornamento del run, per evitare problemi
1891          // dovuti agli eventi annidati (NB! NEVENTS conta anche questi!!)
1892          //
1893          bool fromfirst = true;
1894          //
1895          while ( !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME) && irun < run->GetEntries()-1LL ){
1896    //      while ( iev > (runfirstentry+(ULong64_t)(this->GetRunInfo()->NEVENTS-1)) && irun < run->GetEntries() ){
1897              irun++;
1898              run->GetEntry(irun);
1899              runfirstentry = runlastentry;
1900              if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runfirstentry += 1LL;
1901              runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1902    //        cout << " ))))) UPDATE RUN INFO (((((  @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1903    //        cout << "runfirstentry "<<runfirstentry<<endl;
1904              //      printf(" iev %llu %u %llu \n",iev,this->GetRunInfo()->NEVENTS,(ULong64_t)(runfirstentry+(ULong64_t)(this->GetRunInfo()->NEVENTS)));
1905              //      printf(" abstime %u trailertime %u \n",GetOrbitalInfo()->absTime,GetRunInfo()->RUNTRAILER_TIME);
1906              //      printf(" IDRUN %u \n",GetRunInfo()->ID);
1907              //
1908    //        prevshift = 0;
1909              //
1910              if ( irun == (Long64_t)(run->GetEntries()-1LL) && fromfirst && !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)){
1911                  printf(" resetting irun  (it should NOT happen!!!)\n");
1912                  fromfirst = false;
1913                  irun = 0;
1914                  run->GetEntry(irun);
1915                  runfirstentry = 0ULL;
1916                  runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1917                  if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runlastentry -= 1LL;
1918              };
1919              //
1920          };
1921          //
1922          if (
1923               !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME &&
1924                 GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)
1925              ) {
1926              printf(" Something very wrong here: cannot find RUN containing absolute time %u \n",GetOrbitalInfo()->absTime);
1927              return false;
1928          }
1929          //
1930          if ( irun == oldrun || irun >= run->GetEntries() ) return(false);
1931          //
1932          //  printf(" iev %llu irun %i nevents %u 1st %llu last %llu \n",iev,irun,this->GetRunInfo()->NEVENTS,(ULong64_t)runfirstentry,(ULong64_t)runlastentry);
1933          //
1934          prevshift = 0;
1935          cout << " ))))) UPDATE RUN INFO (((((  @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1936    //      cout << "runfirstentry "<<runfirstentry<<endl;
1937          return(true);    
1938      };
1939      // ----------------------------------------------------
1940      // if it is a preselected file (there is SelectionList)
1941      // NBNB - the event tree MUST be read first
1942      // ----------------------------------------------------
1943      if(SELLI==1){      
1944          sel_tree->GetEntry(iev);
1945    //      cout << irun << " "<< irunentry << endl;
1946          if(irun != oldrun){
1947              run->GetEntry(irun);
1948              cout << " ))))) UPDATE RUN INFO (((((  @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1949              prevshift = 0;
1950              return true;
1951          }
1952          return false;
1953      }
1954    
1955      return false;
1956      //
1957    };
1958    
1959    Bool_t PamLevel2::UpdateRunInfo(Long64_t iev){
1960    
1961        if(!run_tree){
1962            cout << " Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev) -- ERROR -- run tree not loaded"<<endl;
1963            return false;
1964        }
1965        if ( run_tree->GetEntries() <= 0 ) {
1966            cout << " Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev) -- ERROR -- run tree is empty"<<endl;
1967            return(false);
1968        }
1969      
1970        Int_t oldrun = irun; // store current run index
1971    
1972        // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
1973        // if it is a full file (not preselected)
1974        // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
1975        if(SELLI==0){
1976    
1977            // ---------------------------------------------------------------
1978            // the absolute time is necessary to relate the event with the run
1979            // ---------------------------------------------------------------
1980            if( !GetOrbitalInfo() ){
1981                cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- ERROR -- missing OrbitalInfo "<<endl;
1982                return(false);
1983            }
1984          
1985            // -----------------------------------------------------------------------
1986            // the first time the routine is called, set run search from the beginning
1987            // -----------------------------------------------------------------------
1988            if ( irun < 0 ){
1989                irun = 0LL;
1990                run_tree->GetEntry(irun);
1991            };      
1992            //
1993            bool fromfirst = true; // first loop over runs
1994            // ------------------------------------------------------
1995            // loop over runs to find the one that contains the event
1996            // ------------------------------------------------------
1997            while ( !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME &&    // check on absolute time (s)
1998                      GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME) &&
1999                    !(GetOrbitalInfo()->OBT >= GetRunInfo()->RUNHEADER_OBT &&         // additional check on OBT (ms)
2000                      GetOrbitalInfo()->OBT <= GetRunInfo()->RUNTRAILER_OBT) &&
2001                    irun < run_tree->GetEntries() ){
2002                irun++;
2003                // ------------------------------------
2004                // if the end of run tree is reached...
2005                // ------------------------------------
2006                if( irun == run_tree->GetEntries() ){
2007                    if(!fromfirst){
2008                        // -----------------------------------------------------
2009                        // if it happened already once and the run was not found
2010                        // ---> exit with error
2011                        // -----------------------------------------------------
2012                        cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- ERROR -- event entry #"<<iev<<" does not belong to any run (should not happen)" <<endl;
2013                        return false;
2014                    }
2015                    // -----------------------------------------
2016                    // ...otherwise repeat search from beginning
2017                    // -----------------------------------------
2018                    cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- WARNING -- reached end of run tree. searchin again from beginning " <<endl;
2019                    fromfirst = false;
2020                    irun = 0LL;
2021                    runfirstentry = 0LL;
2022                }
2023                // -------------------------------------------------------------------
2024                // save the index of the first entry of the run, relative to pam_tree,
2025                // and read a new run
2026                // -------------------------------------------------------------------
2027                if(irun>0)runfirstentry += (GetRunInfo()->NEVENTS);
2028                run_tree->GetEntry(irun);      
2029                if(GetRunInfo()->RUNHEADER_OBT>=GetRunInfo()->RUNTRAILER_OBT ){
2030                    cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- WARNING -- irun "<<irun<<"  has RUNHEADER_OBT>=RUNTRAILER_OBT " <<endl;
2031                    cout << "                                                            (NB!! in this case some events are assigned to a wrong run)"<<endl;
2032                }
2033            };
2034            // --------------------------------------
2035            // if there was no need to update the run
2036            // ---> exit with FALSE
2037            // --------------------------------------
2038            if ( irun == oldrun ) return(false);
2039    
2040            // --------------------------------------
2041            // ... otherwise
2042            // ---> exit with TRUE
2043            // --------------------------------------
2044            cout << endl << " ))))) UPDATE RUN INFO (((((  @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
2045            prevshift = 0;
2046            return(true);    
2047        };
2048        // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
2049        // if it is a preselected file (there is SelectionList)
2050        // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
2051        if(SELLI==1){      
2052            sel_tree->GetEntry(iev);
2053            if(irun != oldrun){
2054                run_tree->GetEntry(irun);
2055                cout << endl << " ))))) UPDATE RUN INFO (((((  @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
2056                prevshift = 0;
2057                return true;
2058            }
2059            return false;
2060        }
2061    
2062        return false;
2063        //
2064    };
2065    /**
2066     * Update the runinfo informations (to be used to have Run infos event by event basis)
2067     * @param run Pointer to the chain/tree which contains run infos
2068     * @return true if a new run has been read, false if it is still the same run
2069     */
2070    Bool_t PamLevel2::UpdateRunInfo(TTree *run, Long64_t iev){
2071      return(UpdateRunInfo((TChain*)run,iev));
2072    };
2073    
2074    //--------------------------------------
2075    //
2076    //
2077    //--------------------------------------
2078    /**
2079     * Set which trees shoul be analysed
2080     * @param detlist TString containing the sequence of trees required
2081    */
2082    void PamLevel2::SetWhichTrees(TString detlist){
2083            
2084        if(detlist.IsNull() || detlist.Contains("+ALL", TString::kIgnoreCase)){
2085            CAL0 = false;
2086            CAL1 = true;
2087            CAL2 = true;
2088            TRK2 = true;
2089            TRK1 = false;
2090            TRKh = false;
2091            TRK0 = false;
2092            TRG = true;
2093            TOF = true;
2094            TOF0 = false;
2095            S4  = true;
2096            ND  = true;
2097            AC  = true;
2098            ORB = true;
2099        }else if( detlist.Contains("-ALL", TString::kIgnoreCase) ){
2100            CAL0 = false;
2101            CAL1 = false;
2102            CAL2 = false;
2103            TRK2 = false;
2104            TRK1 = false;
2105            TRKh = false;
2106            TRK0 = false;
2107            TRG = false;
2108            TOF = false;
2109            TOF0 = false;
2110            S4  = false;
2111            ND  = false;
2112            AC  = false;
2113            ORB = false;
2114        };
2115        
2116    //  -------------------------------------------------------------------------
2117        if( detlist.Contains("CAL1", TString::kIgnoreCase) ){
2118            if ( detlist.Contains("-CAL1", TString::kIgnoreCase) )CAL1=false;
2119            if ( detlist.Contains("+CAL1", TString::kIgnoreCase) )CAL1=true;
2120        };
2121    
2122        if( detlist.Contains("CAL0", TString::kIgnoreCase) ){
2123            if ( detlist.Contains("-CAL0", TString::kIgnoreCase) )CAL0=false;
2124            if ( detlist.Contains("+CAL0", TString::kIgnoreCase) )CAL0=true;
2125        };
2126            
2127        if( detlist.Contains("CAL2", TString::kIgnoreCase)){
2128            if ( detlist.Contains("-CAL2", TString::kIgnoreCase) )CAL2=false;
2129            if ( detlist.Contains("+CAL2", TString::kIgnoreCase) )CAL2=true;
2130        };
2131            
2132        if( detlist.Contains("+CAL", TString::kIgnoreCase) && !CAL1 && !CAL2 )CAL2=true;
2133        if( detlist.Contains("-CAL", TString::kIgnoreCase) && CAL1 && CAL2 ){
2134            CAL2=false;
2135            CAL1=false;
2136        }
2137    //  -------------------------------------------------------------------------
2138        if( detlist.Contains("TRK0", TString::kIgnoreCase) ){
2139            if ( detlist.Contains("-TRK0", TString::kIgnoreCase) )TRK0=false;
2140            if ( detlist.Contains("+TRK0", TString::kIgnoreCase) )TRK0=true;
2141        };
2142    
2143        if( detlist.Contains("TRK1", TString::kIgnoreCase) ){
2144            if ( detlist.Contains("-TRK1", TString::kIgnoreCase) )TRK1=false;
2145            if ( detlist.Contains("+TRK1", TString::kIgnoreCase) )TRK1=true;
2146        };
2147            
2148        if( detlist.Contains("TRK2", TString::kIgnoreCase)){
2149            if ( detlist.Contains("-TRK2", TString::kIgnoreCase) )TRK2=false;
2150            if ( detlist.Contains("+TRK2", TString::kIgnoreCase) )TRK2=true;
2151        };
2152    
2153        if( detlist.Contains("TRKh", TString::kIgnoreCase)){
2154            if ( detlist.Contains("-TRKh", TString::kIgnoreCase) )TRKh=false;
2155            if ( detlist.Contains("+TRKh", TString::kIgnoreCase) )TRKh=true;
2156        };
2157            
2158        if( detlist.Contains("+TRK", TString::kIgnoreCase) && !TRK1 && !TRK2 && !TRKh )TRK2=true;
2159        if( detlist.Contains("-TRK", TString::kIgnoreCase) && TRK1 && TRK2 && TRKh){
2160            TRK2=false;
2161            TRK1=false;
2162            TRKh=false;
2163        }
2164    //  -------------------------------------------------------------------------
2165        
2166        if( detlist.Contains("-TRG", TString::kIgnoreCase) )TRG = false;
2167        else if( detlist.Contains("+TRG", TString::kIgnoreCase) )TRG = true;
2168        
2169        if( detlist.Contains("-TOF", TString::kIgnoreCase) )TOF = false;
2170        else if( detlist.Contains("+TOF", TString::kIgnoreCase) )TOF = true;
2171    
2172        if( detlist.Contains("-TOF0", TString::kIgnoreCase) )TOF0 = false;
2173        else if( detlist.Contains("+TOF0", TString::kIgnoreCase) )TOF0 = true;
2174        
2175        if( detlist.Contains("-S4",  TString::kIgnoreCase) )S4  = false;
2176        else if( detlist.Contains("+S4",  TString::kIgnoreCase) )S4  = true;
2177        
2178        if( detlist.Contains("-ND",  TString::kIgnoreCase) )ND  = false;
2179        else if( detlist.Contains("+ND",  TString::kIgnoreCase) )ND  = true;
2180        
2181        if( detlist.Contains("-AC",  TString::kIgnoreCase) )AC  = false;
2182        else if( detlist.Contains("+AC",  TString::kIgnoreCase) )AC  = true;
2183        
2184        if( detlist.Contains("-ORB", TString::kIgnoreCase) )ORB = false;
2185        else if( detlist.Contains("+ORB", TString::kIgnoreCase) )ORB = true;
2186    
2187    //     cout<< "Set detector list --> ";
2188    //     if(TRK1)cout<<"TRK1 ";
2189    //     if(TRK2)cout<<"TRK2 ";
2190    //     if(TRKh)cout<<"TRKH ";
2191    //     if(CAL1)cout<<"CAL1 ";
2192    //     if(CAL2)cout<<"CAL2 ";
2193    //     if(TOF)cout<<"TOF ";
2194    //     if(TRG)cout<<"TRG ";
2195    //     if(AC)cout<<"AC ";
2196    //     if(ND)cout<<"ND ";
2197    //     if(S4)cout<<"S4 ";
2198    //     if(ORB)cout<<"ORB ";
2199    //     cout << endl;
2200        
2201    };
2202    
2203    
2204    /**
2205     * Set tree/branch detector flags from the content of a tree
2206     */
2207    void  PamLevel2::GetWhichTrees(TFile* f){
2208        
2209    
2210    
2211        cout << "void  PamLevel2::GetWhichTrees(TFile* f) --- WARNING!! --- ...potrebbe non funzionare "<<endl;
2212        // -----------
2213        // reset flags
2214        // -----------
2215        CAL1   = false;    
2216        CAL2   = false;    
2217        TRK2   = false;    
2218        TRK1   = false;    
2219        TRKh   = false;    
2220        TRG    = false;    
2221        TOF    = false;    
2222        S4     = false;    
2223        ND     = false;    
2224        AC     = false;    
2225        ORB    = false;    
2226        
2227        RUN    = false;
2228            
2229        cout << "Checking file: "<<f->GetName()<<endl;
2230        if( !f || f->IsZombie() ){
2231            cout << "File: "<< f->GetName() <<" Non valid root file"<< endl;
2232            return;
2233        }
2234    
2235        TList *lk = f->GetListOfKeys();
2236        TIter next(lk);
2237        TKey *key =0;
2238    
2239        Int_t nev = 0;
2240    
2241        while( (key = (TKey*)next()) ){
2242            
2243            if( !strcmp(key->GetName(),"Run"        ) )RUN = true;
2244    
2245            //=========================================================    
2246            if( !strcmp(key->GetName(),"Trigger"    ) ){
2247                TRG = true;
2248                Int_t nevt = ((TTree*)f->Get("Trigger"))->GetEntries();
2249                if( nev && nevt!=nev){
2250                    cout << "File: "<< f->GetName() <<" Trigger tree has "<<nevt<<" events instead of "<<nev<< endl;
2251                    TRG = false;
2252                }else nev=nevt;
2253            }
2254            //=========================================================    
2255            if( !strcmp(key->GetName(),"ToF"        ) ){
2256                TOF = true;
2257                Int_t nevt = ((TTree*)f->Get("ToF"))->GetEntries();
2258                if( nev && nevt!=nev){
2259                    cout << "File: "<< f->GetName() <<"     ToF tree has "<<nevt<<" events instead of "<<nev<< endl;
2260                    TOF = false;
2261                }else nev=nevt;
2262            }
2263            //=========================================================  
2264            if( !strcmp(key->GetName(),"S4"         ) ){
2265                S4 = true;
2266                Int_t nevt = ((TTree*)f->Get("S4"))->GetEntries();
2267                if( nev && nevt!=nev){
2268                    cout << "File: "<< f->GetName() <<"      S4 tree has "<<nevt<<" events instead of "<<nev<< endl;
2269                    S4 = false;
2270                }else nev=nevt;
2271            }
2272            //=========================================================  
2273    
2274            if( !strcmp(key->GetName(),"NeutronD"   ) ){
2275                ND = true;
2276                Int_t nevt = ((TTree*)f->Get("NeutronD"))->GetEntries();
2277                if( nev && nevt!=nev){
2278                    cout << "File: "<< f->GetName() <<"NeutronD tree has "<<nevt<<" events instead of "<<nev<< endl;
2279                    ND =false;
2280                }else nev=nevt;
2281            }      
2282            //=========================================================  
2283            if( !strcmp(key->GetName(),"Anticounter") ){
2284                AC = true;
2285                Int_t nevt = ((TTree*)f->Get("Anticounter"))->GetEntries();
2286                if( nev && nevt!=nev){
2287                    cout << "File: "<< f->GetName() <<" Anticounter tree has "<<nevt<<" events instead of "<<nev<< endl;
2288                    AC =false;
2289                }else nev=nevt;
2290            }
2291            //=========================================================  
2292            if( !strcmp(key->GetName(),"OrbitalInfo") ){
2293                ORB = true;
2294                Int_t nevt = ((TTree*)f->Get("OrbitalInfo"))->GetEntries();
2295                if( nev && nevt!=nev){
2296                    cout << "File: "<< f->GetName() <<" OrbitalInfo tree has "<<nevt<<" events instead of "<<nev<< endl;
2297                    ORB = false;
2298                }else nev=nevt;
2299            }
2300            //=========================================================  
2301            if( !strcmp(key->GetName(),"Tracker"    ) ){
2302                TTree *T = (TTree*)f->Get("Tracker");
2303                for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2304                    TString name = T->GetListOfBranches()->At(i)->GetName();
2305                    if( !name.CompareTo("TrkLevel1") )TRK1=true;
2306                    if( !name.CompareTo("TrkLevel2") )TRK2=true;
2307                    if( !name.CompareTo("TrkHough") )TRKh=true;
2308                };      
2309                Int_t nevt = T->GetEntries();
2310                if( nev && nevt!=nev){
2311                    cout << "File: "<< f->GetName() <<" Tracker tree has "<<nevt<<" events instead of "<<nev<< endl;
2312                    TRK1 = false;
2313                    TRK2 = false;
2314                    TRKh = false;
2315                }else nev=nevt;
2316                T->Delete();
2317            };
2318            //=========================================================  
2319            if( !strcmp(key->GetName(),"Calorimeter"    ) ){
2320                TTree *T = (TTree*)f->Get("Calorimeter");
2321                for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2322                    TString name = T->GetListOfBranches()->At(i)->GetName();
2323                    if( !name.CompareTo("CaloLevel1") )CAL1=true;
2324                    if( !name.CompareTo("CaloLevel2") )CAL2=true;
2325                };    
2326                Int_t nevt = T->GetEntries();
2327                if( nev && nevt!=nev){
2328                    cout << "File: "<< f->GetName() <<"  Calorimeter tree has "<<nevt<<" events instead of "<<nev<< endl;
2329                    CAL1 = false;
2330                    CAL2 = false;
2331                }else nev=nevt;
2332                T->Delete();
2333            };      
2334    
2335        };
2336        
2337        delete lk;
2338    
2339    //     cout<< "Get detector list from input file --> ";
2340    //     if(TRK1)cout<<"TRK1 ";
2341    //     if(TRK2)cout<<"TRK2 ";
2342    //     if(TRKh)cout<<"TRKH ";
2343    //     if(CAL1)cout<<"CAL1 ";
2344    //     if(CAL2)cout<<"CAL2 ";
2345    //     if(TOF)cout<<"TOF ";
2346    //     if(TRG)cout<<"TRG ";
2347    //     if(AC)cout<<"AC ";
2348    //     if(ND)cout<<"ND ";
2349    //     if(S4)cout<<"S4 ";
2350    //     if(ORB)cout<<"ORB ";
2351    //     cout << endl;
2352          
2353        return ;
2354            
2355    };
2356    
2357    
2358    //--------------------------------------
2359    //
2360    //
2361    //--------------------------------------
2362    /**
2363     * Check if a file contains selected Pamela Level2 trees.
2364     * @param name File name
2365     * @return true if the file is ok.
2366     */
2367    Bool_t  PamLevel2::CheckLevel2File(TString name){
2368            
2369        Bool_t CAL1__ok   = false;    
2370        Bool_t CAL2__ok   = false;    
2371        Bool_t TRK2__ok   = false;    
2372        Bool_t TRK1__ok   = false;    
2373        Bool_t TRKh__ok   = false;    
2374        Bool_t TRG__ok    = false;    
2375        Bool_t TOF__ok    = false;    
2376        Bool_t S4__ok     = false;    
2377        Bool_t ND__ok     = false;    
2378        Bool_t AC__ok     = false;    
2379        Bool_t ORB__ok    = false;    
2380        
2381        Bool_t RUN__ok    = false;
2382        
2383        Bool_t SELLI__ok  = false;
2384    
2385        cout << "Checking file: "<<name<<endl;
2386        TFile *f = new TFile(name.Data());
2387        if( !f || f->IsZombie() ){
2388            cout << "File: "<< f->GetName() <<" discarded ---- Non valid root file"<< endl; return false;
2389        }
2390    //    cout << "Get list of keys: "<<f<<endl;
2391        TList *lk = f->GetListOfKeys();
2392    //    lk->Print();
2393        TIter next(lk);
2394        TKey *key =0;
2395    
2396        Int_t nev = 0;
2397    
2398        while( (key = (TKey*)next()) ){
2399            
2400    //      cout << key->GetName() << endl;
2401    //      cout << key->GetName() << ""<<key->GetClassName()<<endl;
2402    //              cout << " Get tree: " << f->Get(key->GetName())<<endl;
2403    //      nev_previous = nev;
2404    //      cout << " n.entries  "<< nev <<endl;
2405    //      if( key->GetClassName()=="TTree" && nev_previous && nev != nev_previous ){
2406    //          nev = ((TTree*)f->Get(key->GetName()))->GetEntries();
2407    //          cout << "File: "<< f->GetName() <<" discarded ---- "<< key->GetName() << " tree: n.entries does not match "<<nev<<" "<<nev_previous<< endl;
2408    //          return false;
2409    //      };
2410    
2411    
2412            if( !strcmp(key->GetName(),"Run"        ) )RUN__ok = true;      
2413    
2414            //=========================================================    
2415            if( !strcmp(key->GetName(),"SelectionList"    ) ){
2416                SELLI__ok = true;
2417                if(SELLI==1){
2418                    Int_t nevt = ((TTree*)f->Get("SelectionList"))->GetEntries();
2419                    if( nev && nevt!=nev){
2420                        cout << "File: "<< f->GetName() <<" discarded ---- SelectionList tree has "<<nevt<<" events instead of "<<nev<< endl;
2421                        return false;
2422                    }
2423                    nev=nevt;
2424                }
2425            }
2426    
2427            //=========================================================    
2428            if( !strcmp(key->GetName(),"Trigger"    ) ){
2429                TRG__ok = true;
2430                if(TRG){
2431                    Int_t nevt = ((TTree*)f->Get("Trigger"))->GetEntries();
2432                    if( nev && nevt!=nev){
2433                        cout << "File: "<< f->GetName() <<" discarded ---- Trigger tree has "<<nevt<<" events instead of "<<nev<< endl;
2434                        return false;
2435                    }
2436                    nev=nevt;
2437                }
2438            }
2439            //=========================================================    
2440            if( !strcmp(key->GetName(),"ToF"        ) ){
2441                TOF__ok = true;
2442                if(TOF){
2443                    Int_t nevt = ((TTree*)f->Get("ToF"))->GetEntries();
2444                    if( nev && nevt!=nev){
2445                        cout << "File: "<< f->GetName() <<" discarded ---- ToF tree has "<<nevt<<" events instead of "<<nev<< endl;
2446                        return false;
2447                    }
2448                    nev=nevt;
2449                }
2450            }
2451            //=========================================================  
2452            if( !strcmp(key->GetName(),"S4"         ) ){
2453                S4__ok = true;
2454                if(S4){
2455                    Int_t nevt = ((TTree*)f->Get("S4"))->GetEntries();
2456                    if( nev && nevt!=nev){
2457                        cout << "File: "<< f->GetName() <<" discarded ---- S4 tree has "<<nevt<<" events instead of "<<nev<< endl;
2458                        return false;
2459                    }
2460                    nev=nevt;
2461                }
2462            }
2463            //=========================================================  
2464    
2465            if( !strcmp(key->GetName(),"NeutronD"   ) ){
2466                ND__ok = true;
2467                if(ND){
2468                    Int_t nevt = ((TTree*)f->Get("NeutronD"))->GetEntries();
2469                    if( nev && nevt!=nev){
2470                        cout << "File: "<< f->GetName() <<" discarded ---- NeutronD tree has "<<nevt<<" events instead of "<<nev<< endl;
2471                        return false;
2472                    }
2473                    nev=nevt;
2474                }
2475            }      
2476            //=========================================================  
2477            if( !strcmp(key->GetName(),"Anticounter") ){
2478                AC__ok = true;
2479                if(AC){
2480                    Int_t nevt = ((TTree*)f->Get("Anticounter"))->GetEntries();
2481                    if( nev && nevt!=nev){
2482                        cout << "File: "<< f->GetName() <<" discarded ---- Anticounter tree has "<<nevt<<" events instead of "<<nev<< endl;
2483                        return false;
2484                    }
2485                    nev=nevt;
2486                }
2487            }
2488            //=========================================================  
2489            if( !strcmp(key->GetName(),"OrbitalInfo") ){
2490                ORB__ok = true;
2491                if(ORB){
2492                    Int_t nevt = ((TTree*)f->Get("OrbitalInfo"))->GetEntries();
2493                    if( nev && nevt!=nev){
2494                        cout << "File: "<< f->GetName() <<" discarded ---- OrbitalInfo tree has "<<nevt<<" events instead of "<<nev<< endl;
2495                        return false;
2496                    }
2497                    nev=nevt;
2498                }
2499            }
2500            //=========================================================  
2501            if( !strcmp(key->GetName(),"Tracker"    ) ){
2502                TTree *T = (TTree*)f->Get("Tracker");
2503                if(TRK1||TRK2||TRKh){
2504                    Int_t nevt = T->GetEntries();
2505                    if( nev && nevt!=nev){
2506                        cout << "File: "<< f->GetName() <<" discarded ---- Tracker tree has "<<nevt<<" events instead of "<<nev<< endl;
2507                        return false;
2508                    }
2509                    nev=nevt;
2510                }
2511                for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2512                    TString name = T->GetListOfBranches()->At(i)->GetName();
2513                    if( !name.CompareTo("TrkLevel1") )TRK1__ok=true;
2514                    if( !name.CompareTo("TrkLevel2") )TRK2__ok=true;
2515                    if( !name.CompareTo("TrkHough") )TRKh__ok=true;
2516                };      
2517                T->Delete();
2518            };
2519            //=========================================================  
2520            if( !strcmp(key->GetName(),"Calorimeter"    ) ){
2521                TTree *T = (TTree*)f->Get("Calorimeter");
2522                if(CAL1||CAL2){
2523                    Int_t nevt = T->GetEntries();
2524                    if( nev && nevt!=nev){
2525                        cout << "File: "<< f->GetName() <<" discarded ---- Calorimeter tree has "<<nevt<<" events instead of "<<nev<< endl;
2526                        return false;
2527                    }
2528                    nev=nevt;
2529                }
2530                for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2531                    TString name = T->GetListOfBranches()->At(i)->GetName();
2532                    if( !name.CompareTo("CaloLevel1") )CAL1__ok=true;
2533                    if( !name.CompareTo("CaloLevel2") )CAL2__ok=true;
2534                };    
2535                T->Delete();
2536            };      
2537    
2538        };
2539    
2540        if( SELLI==-1 )SELLI = (Int_t)SELLI__ok;
2541        if( SELLI==0 && SELLI__ok ){
2542            cout << "File: "<< f->GetName() <<" discarded ---- found SelectionList (it is not a full-event file)" << endl;
2543            return false;  
2544        }
2545        if( SELLI==1 && !SELLI__ok ){
2546            cout << "File: "<< f->GetName() <<" discarded ---- SelectionList missing" << endl;
2547            return false;
2548        }
2549        
2550    //    cout << "SELLI "<<SELLI<<endl;
2551    
2552    //     cout<< "CheckLevel2File(TString): detector list --> ";
2553    //     if(TRK1__ok)cout<<"TRK1 ";
2554    //     if(TRK2__ok)cout<<"TRK2 ";
2555    //     if(TRKh__ok)cout<<"TRKH ";
2556    //     if(CAL1__ok)cout<<"CAL1 ";
2557    //     if(CAL2__ok)cout<<"CAL2 ";
2558    //     if(TOF__ok)cout<<"TOF ";
2559    //     if(TRG__ok)cout<<"TRG ";
2560    //     if(AC__ok)cout<<"AC ";
2561    //     if(ND__ok)cout<<"ND ";
2562    //     if(S4__ok)cout<<"S4 ";
2563    //     if(ORB__ok)cout<<"ORB ";
2564    //     cout << endl;
2565    
2566    
2567        if(TRK2 && TRK1__ok)TRK1=1;
2568    // ----------------------------------------------------------------------------
2569    // NOTA
2570    // se c'e` il level1, lo devo necessarimente leggere.
2571    // infatti (non ho capito perche`) i cluster vengono letti e allocati in memoria
2572    // comunque, ma non vengono disallocati da PamLevel2::Clear()
2573    // ----------------------------------------------------------------------------
2574    
2575    
2576        if(!RUN__ok) {
2577            cout << "File: "<< f->GetName() <<" *WARNING* ---- Missing RunInfo tree"<< endl;
2578    //      return false;
2579        };
2580    
2581        if(CAL1 && !CAL1__ok){
2582            cout << "File: "<< f->GetName() <<" discarded ---- Missing CaloLevel1 branch"<< endl;
2583            return false;
2584        };
2585        if(CAL2 && !CAL2__ok){
2586            cout << "File: "<< f->GetName() <<" discarded ---- Missing CaloLevel2 branch"<< endl;
2587            return false;
2588        };
2589        if(TRK2 && !TRK2__ok){
2590            cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkLevel2 branch"<< endl;
2591            return false;
2592        };
2593        if(TRK1 && !TRK1__ok){
2594            cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkLevel1 branch"<< endl;
2595            return false;
2596        };
2597        if(TRKh && !TRKh__ok){
2598            cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkHough branch"<< endl;
2599            return false;
2600        };
2601        if(ORB && !ORB__ok){
2602            cout << "File: "<< f->GetName() <<" discarded ---- Missing ORB tree"<< endl;
2603            return false;
2604        };
2605        if(AC && !AC__ok){
2606            cout << "File: "<< f->GetName() <<" discarded ---- Missing AC tree"<< endl;
2607            return false;
2608        };
2609        if(S4 && !S4__ok){
2610            cout << "File: "<< f->GetName() <<" discarded ---- Missing S4 tree"<< endl;
2611            return false;
2612        };
2613        if(TOF && !TOF__ok){
2614            cout << "File: "<< f->GetName() <<" discarded ---- Missing ToF tree"<< endl;
2615            return false;
2616        };
2617    
2618        if(ND && !ND__ok){
2619            cout << "File: "<< f->GetName() <<" discarded ---- Missing ND tree"<< endl;
2620            return false;
2621        };
2622        if(TRG && !TRG__ok){
2623            cout << "File: "<< f->GetName() <<" discarded ---- Missing Trigger tree"<< endl;
2624            return false;
2625        };
2626        
2627    
2628    //    lk->Delete();
2629    //    delete lk;
2630        f->Close();
2631    
2632    //     cout<< "CheckLevel2File(TString): detector list --> ";
2633    //     if(TRK1)cout<<"TRK1 ";
2634    //     if(TRK2)cout<<"TRK2 ";
2635    //     if(TRKh)cout<<"TRKH ";
2636    //     if(CAL1)cout<<"CAL1 ";
2637    //     if(CAL2)cout<<"CAL2 ";
2638    //     if(TOF)cout<<"TOF ";
2639    //     if(TRG)cout<<"TRG ";
2640    //     if(AC)cout<<"AC ";
2641    //     if(ND)cout<<"ND ";
2642    //     if(S4)cout<<"S4 ";
2643    //     if(ORB)cout<<"ORB ";
2644    //     cout << endl;
2645          
2646        return true;
2647            
2648    };
2649    
2650    
2651    /**
2652     * Create clone-trees
2653     */
2654    void PamLevel2::CreateCloneTrees0( TChain *fChain, TFile *ofile ){
2655    
2656        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2657        cout << "Create clones of PAMELA trees "<<endl;
2658        
2659        Int_t i=0;
2660        pam_tree_clone[i] = fChain->GetTree()->CloneTree(0);
2661        TString name =  pam_tree_clone[i]->GetName();
2662        name.Append("_clone");
2663    //    pam_tree_clone[i]->SetName(name.Data());
2664        cout << pam_tree_clone[i]->GetName() <<endl;
2665        i++;
2666    
2667        TList *li = fChain->GetListOfFriends();
2668        TIter next(li);
2669        TFriendElement* T_friend=0;
2670        ofile->cd();
2671        while( (T_friend = (TFriendElement*)next()) ){
2672    //      cout<<T_friend->IsA()->GetName()<<" "<<T_friend->GetName()<<hex << T_friend->GetTree() << dec<<endl;
2673    //      cout<<T_friend->GetTree()->GetName()<< endl;
2674            pam_tree_clone[i] = T_friend->GetTree()->CloneTree(0);
2675            pam_tree_clone[i]->SetAutoSave(1000000);
2676            name =  pam_tree_clone[i]->GetName();
2677            name.Append("_clone");
2678    //      pam_tree_clone[i]->SetName(name.Data());
2679            cout << pam_tree_clone[i]->GetName() << endl;
2680            i++;
2681        }
2682        
2683        delete li;
2684    
2685        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2686    
2687    }
2688    
2689    /**
2690     * Create clone-trees
2691     */
2692    void PamLevel2::CreateCloneTrees(TFile *ofile){
2693    
2694        ofile->cd();
2695    
2696        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2697        cout << "Create new PAMELA trees "<<endl;
2698            
2699    
2700        run_tree_clone = new TTree("Run","PAMELA Level2 data from the GL_RUN table ");
2701        run_tree_clone->Branch("RunInfo","GL_RUN",GetPointerTo("RunInfo"));
2702        cout << "Run          : branch RunInfo"<<endl;
2703        run_tree_clone->Branch("SoftInfo","SoftInfo",GetPointerTo("SoftInfo"));
2704        cout << "Run          : branch SoftInfo"<<endl;
2705        // ------------------
2706        // replicate run tree
2707        // ------------------
2708        cout << "----------------------------------------------------"<<endl;
2709        cout << "irun\t | RUN\t NEVENTS\t absolute time"<<endl;
2710        for (Int_t i=0; i<run_tree->GetEntries(); i++){
2711            run_tree->GetEntry(i);
2712            cout << i<< "\t | "<<GetRunInfo()->ID<<"\t "<<GetRunInfo()->NEVENTS<< "\t "<<GetRunInfo()->RUNHEADER_TIME<<" <---> "<<GetRunInfo()->RUNTRAILER_TIME<<endl;
2713            run_tree_clone->Fill();
2714        }
2715        cout << "----------------------------------------------------"<<endl;
2716    
2717    
2718        sel_tree_clone = new TTree("SelectionList","List of selected events ");
2719        sel_tree_clone->Branch("RunEntry",&irun,"runentry/L");
2720        sel_tree_clone->Branch("EventEntry",&irunentry,"eventry/L");
2721        
2722    
2723        Int_t i=0;
2724        if(TRK1||TRK2||TRKh){
2725            pam_tree_clone[i] = new TTree("Tracker","PAMELA tracker level2 data ");
2726            if(TRK1) {
2727                pam_tree_clone[i]->Branch("TrkLevel1","TrkLevel1", GetPointerTo("TrkLevel1"));
2728                pam_tree_clone[i]->BranchRef();
2729                cout << "Tracker      : branch TrkLevel1"<<endl;
2730    //          cout << "CreateCloneTrees " << GetTrkLevel1()<<endl;
2731            };
2732            if(TRK2) {
2733                pam_tree_clone[i]->Branch("TrkLevel2", "TrkLevel2",GetPointerTo("TrkLevel2"));
2734                cout << "Tracker      : branch TrkLevel2"<<endl;        
2735            };
2736            if(TRKh) {
2737                pam_tree_clone[i]->Branch("TrkHough","TrkHough", GetPointerTo("TrkHough"));
2738                cout << "Tracker      : branch TrkHough"<<endl;
2739            };
2740            i++;
2741        }
2742    
2743        // Calorimeter
2744        if(CAL1||CAL2){
2745            pam_tree_clone[i] = new TTree("Calorimeter","PAMELA calorimeter level2 data ");
2746            if(CAL1) {
2747                pam_tree_clone[i]->Branch("CaloLevel1", "CaloLevel1", GetPointerTo("CaloLevel1"));
2748                cout << "Calorimeter  : branch CaloLevel1"<<endl;
2749            };
2750            if(CAL2) {
2751                pam_tree_clone[i]->Branch("CaloLevel2","CaloLevel2", GetPointerTo("CaloLevel2"));
2752                cout << "Calorimeter  : branch CaloLevel2"<<endl;
2753            };
2754            i++;
2755        }    
2756    
2757        // ToF    
2758        if(TOF) {
2759            pam_tree_clone[i] = new TTree("ToF","PAMELA ToF level2 data ");
2760            pam_tree_clone[i]->Branch("ToFLevel2","ToFLevel2", GetPointerTo("ToFLevel2"));
2761            cout << "ToF          : branch ToFLevel2"<<endl;
2762            i++;
2763        };
2764        // Trigger
2765        if(TRG) {
2766            pam_tree_clone[i] = new TTree("Trigger","PAMELA trigger level2 data ");
2767            pam_tree_clone[i]->Branch("TrigLevel2","TrigLevel2", GetPointerTo("TrigLevel2"));
2768            cout << "Trigger      : branch TrigLevel2"<<endl;
2769            i++;
2770        };
2771        // S4
2772        if(S4) {
2773            pam_tree_clone[i] = new TTree("S4","PAMELA S4 level2 data ");  
2774            pam_tree_clone[i]->Branch("S4Level2","S4Level2", GetPointerTo("S4Level2"));
2775            cout << "S4           : branch S4Level2"<<endl;
2776            i++;
2777        };
2778        // Neutron Detector
2779        if(ND) {
2780            pam_tree_clone[i] = new TTree("NeutronD","PAMELA neutron detector level2 data ");      
2781            pam_tree_clone[i]->Branch("NDLevel2","NDLevel2", GetPointerTo("NDLevel2"));
2782            cout << "NeutronD     : branch NDLevel2"<<endl;
2783            i++;
2784        };
2785        // Anticounters
2786        if(AC) {
2787            pam_tree_clone[i] = new TTree("Anticounter","PAMELA anticounter detector level2 data ");        
2788            pam_tree_clone[i]->Branch("AcLevel2","AcLevel2", GetPointerTo("AcLevel2"));
2789            cout << "Anticounter  : branch AcLevel2"<<endl;
2790            i++;
2791        };
2792        // OrbitalInfo
2793        if(ORB) {
2794            pam_tree_clone[i] = new TTree("OrbitalInfo","PAMELA oribital info  ");  
2795            pam_tree_clone[i]->Branch("OrbitalInfo","OrbitalInfo", GetPointerTo("OrbitalInfo"));
2796            cout << "OrbitalInfo  : branch OrbitalInfo"<<endl;
2797            i++;
2798        };
2799        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2800    
2801    }
2802    
2803    
2804    /**
2805     * Fill tree (created with CreateCloneTrees)
2806     *
2807     */
2808    //void PamLevel2::FillNewPamTree(TTree *T){
2809    void PamLevel2::FillCloneTrees(){
2810        
2811    //    cout << "PamLevel2::FillCloneTrees()" << irunentry << endl;
2812    
2813        for(Int_t i=0; i<NCLONES; i++){
2814            if(pam_tree_clone[i])pam_tree_clone[i]->Fill();
2815        }
2816        if(sel_tree_clone)sel_tree_clone->Fill();    
2817    
2818    }
2819    
2820    
2821    TTree* PamLevel2::GetCloneTree(TString name){
2822    
2823        for(Int_t i=0; i<NCLONES; i++){    
2824            if(pam_tree_clone[i]){
2825                TString na = pam_tree_clone[i]->GetName();
2826                if(!name.CompareTo(na))return pam_tree_clone[i];
2827            };
2828        }
2829        if(run_tree_clone){
2830            TString na = run_tree_clone->GetName();
2831            if(!name.CompareTo(na))return run_tree_clone;
2832        }
2833        if(sel_tree_clone){
2834            TString na = sel_tree_clone->GetName();
2835            if(!name.CompareTo(na))return sel_tree_clone;
2836        }
2837        return NULL;
2838    
2839    }
2840    void PamLevel2::WriteCloneTrees(){
2841        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2842        cout << "Write clones of PAMELA trees "<<endl;
2843        cout << run_tree_clone->GetName()<<endl;    
2844        run_tree_clone->Write();
2845        cout << sel_tree_clone->GetName()<<endl;    
2846        sel_tree_clone->Write();
2847        for(Int_t i=0; i<NCLONES; i++){    
2848            if(pam_tree_clone[i]){
2849                cout << pam_tree_clone[i]->GetName()<<endl;
2850                pam_tree_clone[i]->Write();
2851            };
2852        }
2853        cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2854    
2855    }
2856    
2857    /**
2858     * Method to get level2-trees entry, the corresponding run entry and (if required) the level0 entry.
2859     */
2860    //Int_t PamLevel2::GetEntry(Int_t iee){
2861    Int_t PamLevel2::GetEntry(Long64_t iee){
2862        
2863        if(!pam_tree){
2864            cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- level2 trees not loaded"<<endl;
2865            return 0;
2866        }
2867    
2868    
2869        //
2870        // This is a sort of bug: if you don't have the run tree you don't want to exit here you want to have loaded the event anyway...
2871        //
2872        //    if(!run_tree ){
2873        //  cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- run tree not loeaded"<<endl;
2874        //  return 0;
2875        //    }
2876    
2877        Long64_t ii=0;
2878        //-------------------------------
2879        ii = iee;
2880        if( !pam_tree->GetEntry(ii) ){      
2881            cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- error reading pam tree"<<endl;
2882            return 0;
2883        }
2884        //
2885        // ... that's way I put it here. Notice that nothing change in the code (is backward compatible) since in any case you return with 0.
2886        // in theory one would like to return 1 if run is not loaded but now I don't have the will to add that 2 lines of code and it is not
2887        // a problem if you don't check the return code of getentry.
2888        //
2889        if(!run_tree ){
2890            if ( TRK0 || CAL0 || TOF0 || RUN ) { //forse cosi` va bene per tornare 1?
2891                cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- run tree not loaded"<<endl;
2892                return 0;
2893            }  else {
2894                return 1; //cosi` se non c'e` run esce qua...
2895            }
2896        }
2897    
2898        //-------------------------------
2899        ii = iee;
2900    //    Bool_t UPDATED = UpdateRunInfo(run_tree,ii);
2901    //    Bool_t UPDATED = UpdateRunInfo(ii);
2902        UpdateRunInfo(ii);
2903        if(SELLI==0)irunentry = iee-runfirstentry;
2904    //    if(UPDATED && run_tree_clone)run_tree_clone->Fill();
2905      
2906    //    cout << "PamLevel2::GetEntry("<<iee<<") "<<irun<<" "<<runfirstentry<<" "<<irunentry<<endl;
2907    
2908    //    cout << " irunentry "<<irunentry << endl;
2909    
2910        if( TRK0 || CAL0 || TOF0 ){
2911            if( !GetYodaEntry( ) ){
2912                cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- error reading level0 tree"<<endl;
2913                return 0;
2914            }
2915        }
2916        return 1;
2917    
2918    }
2919    
2920    /**
2921     * Method to retrieve the level0 tree (YODA tree) that contains the current event.
2922     * Given the run ID (...), if needed it query the DB and load the proper file.
2923     * @return Pointer to the tree
2924     */
2925    
2926    
2927    TTree* PamLevel2::GetYodaTree( ){
2928        
2929    //    cout << "TTree* PamLevel2::GetYodaTree( )"<<endl;
2930        //===================================
2931        // check if iroot has changed
2932        //===================================
2933        if( irun<0 ){
2934            cout << "PamLevel2::GetYodaTree() -- ERROR -- irun = "<<irun<<endl;
2935    //      cout << "In order to use this method you have to first load the RunInfo tree "<<endl;
2936            return NULL;
2937        }    
2938        Int_t irootnew = run_obj->ID_ROOT_L0;
2939    //     cout << "iroot    "<<iroot<<endl;
2940    //     cout << "irootnew "<<irootnew<<endl;
2941    
2942        //===================================
2943        // load the level0 file
2944        // (if not already loaded)
2945        //===================================
2946        if( iroot != irootnew || !l0_tree){
2947            iroot = irootnew;
2948            //===================================
2949            // open the DB connection
2950            // (if not already opened)
2951            //===================================
2952            if(!dbc || (dbc && !dbc->IsConnected())){
2953                cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
2954                cout<<"Connecting to DB"<<endl;
2955                cout<<"HOST "<<host<<endl;
2956                cout<<"USER "<<user<<endl;
2957                cout<<"PSW  "<<psw<<endl;
2958                dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
2959                if( !dbc )return NULL;
2960                if( !dbc->IsConnected() )return NULL;    
2961    //          cout<<"...done"<<endl;
2962                cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
2963            }else{
2964    //          cout<<"DB already connected"<<endl;
2965            }
2966            GL_ROOT glroot = GL_ROOT();
2967            if( glroot.Query_GL_ROOT(iroot,dbc) ){
2968                cout << "TTree* PamLevel2::GetYodaTree( ) -- ERROR -- level0 file iroot = "<<iroot<< " does not exists"<<endl;
2969                return NULL;
2970            };
2971            TString filename = glroot.PATH + glroot.NAME;
2972            if(l0_file){
2973                l0_file->Close();
2974                l0_file->Delete();
2975            }
2976            cout << "Opening LEVEL0 file: "<< filename << endl;
2977            FileStat_t t;
2978            if( gSystem->GetPathInfo(filename.Data(),t) ){
2979                cout << " PamLevel2::GetYodaTree() -- ERROR opening file "<<endl;
2980                return NULL;
2981            }
2982            l0_file = new TFile(filename);
2983            if( !l0_file )return NULL;
2984            l0_tree = (TTree*)l0_file->Get("Physics");
2985            if(!h0_obj)h0_obj = new EventHeader();
2986            l0_tree->SetBranchAddress("Header" ,&h0_obj);
2987            prevshift = 0;
2988            //---------------------------------------------------
2989            // TRACKER:
2990            if(TRK0){
2991                if(!trk0_obj){
2992                    trk0_obj = new TrkLevel0();
2993                    trk0_obj->Set();
2994                };
2995                l0_tree->SetBranchAddress("Tracker" ,trk0_obj->GetPointerToTrackerEvent());
2996            }
2997            //---------------------------------------------------
2998            // CALORIMETER:
2999            if(CAL0){
3000                if(!calo0_obj){
3001                    calo0_obj = new CaloLevel0();
3002                    calo0_obj->Set();
3003                };
3004                l0_tree->SetBranchAddress("Calorimeter" ,calo0_obj->GetPointerToCalorimeterEvent());
3005                //      cout << "PamLevel2::GetYodaTree() --- level0 calorimeter not implemented "<<endl;
3006            }
3007            //---------------------------------------------------
3008            // TOF:
3009            if(TOF0){
3010                cout << "PamLevel2::GetYodaTree() --- level0 TOF not implemented "<<endl;
3011            }
3012    
3013        };
3014    
3015        if(!dbc || (dbc && !dbc->IsConnected())){
3016            cout << " TTree* PamLevel2::GetYodaTree( ) -- no DB connected... hai fatto qualche cazzata "<<endl;
3017        }
3018    
3019        if ( TRK0 ){
3020            TrkParams::Load(6);
3021            if( !TrkParams::IsLoaded(6) ){
3022                cout << " TTree* PamLevel2::GetYodaTree( ) -- WARNING -- VK-mask not loaded"<<endl;
3023            };
3024            TrkParams::SetCalib(run_obj,dbc);
3025            TrkParams::LoadCalib( );
3026            if( !TrkParams::CalibIsLoaded() ){
3027                cout << " TTree* PamLevel2::GetYodaTree( ) -- WARNING -- Calibration not loaded"<<endl;
3028            };
3029        }
3030    
3031    //    cout << l0_tree << endl;
3032        
3033        return l0_tree;
3034    
3035    }
3036    
3037    /**
3038     * Method to retrieve the level0 tree (YODA tree) that contains the current event.
3039     */
3040    Int_t PamLevel2::GetYodaEntry(){
3041    
3042    //    cout << "Int_t PamLevel2::GetYodaEntry()"<<endl;
3043        if(!GetYodaTree())return 0;
3044        
3045        // patch
3046        if( irunentry < 0){
3047    //      cout << "Int_t PamLevel2::GetYodaEntry() -- ATTENZIONE -- irunentry negativo?!?! "<<(Int_t)irunentry<<endl;
3048            irunentry=0LL;
3049        }
3050        //  ---------------------------------
3051        //  if file is NOT a preselected file
3052        //  ---------------------------------
3053        Long64_t quellagiusta = irunentry + (Long64_t)(run_obj->EV_FROM);
3054    //     cout << " irunentry "<<irunentry<<endl;
3055    //     cout << " EV_FROM "<<run_obj->EV_FROM<<endl;
3056    //     cout << " quellagiusta = irunentry + EV_FROM "<< quellagiusta << endl;
3057    
3058    //     cout << " irun "<< irun << " irunentry "<< irunentry<<" run_obj->EV_FROM "<<run_obj->EV_FROM <<" quella giusta "<<quellagiusta << endl;
3059    //     cout << " iroot "<<iroot<<" run_obj->ID_ROOT_L0 "<<run_obj->ID_ROOT_L0<<endl;
3060    //     cout << " time "<< GetOrbitalInfo()->absTime << endl;
3061    //     cout << " trk_calib_used "<<run_obj->TRK_CALIB_USED<< endl;
3062        
3063        if( !GetOrbitalInfo() ){
3064            cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- missing OrbitalInfo "<<endl;
3065            return 0;
3066        }
3067        if( GetOrbitalInfo()->OBT==0 && GetOrbitalInfo()->pkt_num==0 ){
3068            cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- level2 event corrupted ?? "<<endl;
3069            return 0;
3070        }
3071    
3072        // ---------------------------------------------------------------------
3073        // ATTENTION!!!
3074        // If data are not pre-processed with cleaner, the level0 tree may contain
3075        // spurious nested physics packets.
3076        // The GL_RUN variables EV_FROM, EV_TO and NEVENTS counts also these entries
3077        // while level2 tree DOES NOT!!
3078        // This means that "quellagiusta" in these cases is not correct.
3079        // In order to retrieve the correct level0 event, I implemented a check
3080        // of the OBT and pkt-number. In case of mismatch, the level0 entry number
3081        // is shift forward until when the packets match.
3082        // ---------------------------------------------------------------------
3083        Int_t answer = 0;
3084        Int_t shift =0;
3085        //    printf(" siamo qui %i %i \n",shift,prevshift);
3086        Int_t maxshift = 100;
3087        do{
3088            if(shift>0){    
3089                cout << " PKTNUM  L2 --- "<< GetOrbitalInfo()->pkt_num << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetCounter()<<endl;
3090                cout << "         RUN: ID "<< GetRunInfo()->ID << " ID_ROOT_L0 "<<run_obj->ID_ROOT_L0<<" ID_RUN_FRAG "<<GetRunInfo()->ID_RUN_FRAG << " EV_FROM "<<GetRunInfo()->EV_FROM  <<endl;
3091                cout << "         L2 <--> L0 mismatch ( irun "<<irun<<" irunentry "<<irunentry<<" shift "<<shift<<" prevshift "<<prevshift<<" )"<<endl;
3092            }
3093            answer = l0_tree->GetEntry(quellagiusta+(Long64_t)shift+(Long64_t)prevshift);
3094            shift++;
3095            if( !GetEventHeader() ){
3096                cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- missing EventHeader "<<endl;
3097                return 0;
3098            }
3099    //      cout << "PKTNUM "<<shift<<" ==  L2 --- "<< GetOrbitalInfo()->pkt_num << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetCounter()<<endl;
3100    //      cout << " L2 --- "<< GetOrbitalInfo()->OBT << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetOrbitalTime()<<endl;
3101    //      if( (quellagiusta+shift) == l0_tree->GetEntries()+1 )cout << ">>> end of level0 tree <<<"<<endl;
3102    //      cout << " GetOrbitalInfo()->OBT "<< GetOrbitalInfo()->OBT << endl;
3103    //      cout << " GetEventHeader()->GetPscuHeader()->GetOrbitalTime() "<< GetEventHeader()->GetPscuHeader()->GetOrbitalTime() << endl;
3104    //      cout << " GetOrbitalInfo()->pkt_num "<< GetOrbitalInfo()->pkt_num << endl;
3105    //      cout << " GetEventHeader()->GetPscuHeader()->GetCounter() "<< GetEventHeader()->GetPscuHeader()->GetCounter() << endl;
3106    //      printf(" IDRUN %u \n",GetRunInfo()->ID);
3107    //
3108            if ( prevshift != 0 && (quellagiusta+(Long64_t)shift) == GetYodaTree()->GetEntries() ){
3109              prevshift = 0;
3110              shift = -1;
3111            };
3112    
3113        }while( ( GetOrbitalInfo()->OBT != (UInt_t)(GetEventHeader()->GetPscuHeader()->GetOrbitalTime()) || GetOrbitalInfo()->pkt_num != (UInt_t)(GetEventHeader()->GetPscuHeader()->GetCounter())) && (quellagiusta+(Long64_t)shift) < GetYodaTree()->GetEntries() && shift < maxshift);
3114    
3115        if ( (quellagiusta+(Long64_t)shift) == GetYodaTree()->GetEntries() || shift == maxshift ) {
3116            cout << " Big trouble here, no such event in Level0 data! " <<endl;
3117            return 0;
3118        }
3119    //    cout << "LA ENTRY GIUSTA E`: "<<quellagiusta<<" (spero...)"<<endl;
3120    //    return GetYodaTree()->GetEntry(quellagiusta);
3121        if ( shift > 1 ) prevshift += (shift-1);
3122        
3123        return answer;
3124        
3125    }
3126    
3127    

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.37

  ViewVC Help
Powered by ViewVC 1.1.23