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

Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.45

  ViewVC Help
Powered by ViewVC 1.1.23