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

Diff of /PamelaLevel2/src/PamLevel2.cpp

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.23